{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "affecting-percentage",
   "metadata": {},
   "outputs": [],
   "source": [
    "# External libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import os\n",
    "import requests\n",
    "\n",
    "# helpful check routine\n",
    "def is_downloadable(url):\n",
    "    \"\"\"\n",
    "    Does the url contain a downloadable resource\n",
    "    \"\"\"\n",
    "    h = requests.head(url, allow_redirects=True)\n",
    "    header = h.headers\n",
    "    content_type = header.get('content-type')\n",
    "    if 'text' in content_type.lower():\n",
    "        return False\n",
    "    if 'html' in content_type.lower():\n",
    "        return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "asian-frame",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CSI routines\n",
    "import csi.planarfault as pf\n",
    "import csi.gps as gr\n",
    "import csi.insar as insar\n",
    "import csi.geodeticplot as geoplt\n",
    "import csi.seismiclocations as sl\n",
    "import csi.imagedownsampling as down\n",
    "import csi.fault3D as flt3D\n",
    "import csi.faultwithvaryingdip as flt\n",
    "import csi.verticalfault as verticalfault\n",
    "#import csi.cosicorrrates as cr\n",
    "import csi.imagecovariance as imcov\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "received-spiritual",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set some flags\n",
    "input_check = True\n",
    "do_downsample = False\n",
    "output_check = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "developed-fence",
   "metadata": {},
   "source": [
    "CSI does most computations in a local Cartesian coordinate system. It presently supports only UTM projections, so we need to set the zone to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "abstract-tanzania",
   "metadata": {},
   "outputs": [],
   "source": [
    "# UTM zone 11 for eastern California\n",
    "utmzone = 11"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fatal-japan",
   "metadata": {},
   "source": [
    "For this example, we will use an InSAR dataset for the 2019 Ridgecrest earthquake, published online at the Harvard Dataverse:\n",
    "\n",
    "Fielding, Eric Jameson, 2019, \"Replication Data for: Surface deformation related to the 2019 Mw 7.1 and Mw 6.4 Ridgecrest Earthquakes in California from GPS, SAR interferometry, and SAR pixel offsets\", https://doi.org/10.7910/DVN/JL9YMS, Harvard Dataverse, V2\n",
    "\n",
    "The dataset we will start with is the ALOS-2 interferogram from ascending track 65. We set the output name for the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "stuffed-capital",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the output name root\n",
    "outName = 'A2A065'  # ALOS-2 Ascending track 65"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prime-momentum",
   "metadata": {},
   "source": [
    "We use the masked range-change file to avoid using data from areas with low coherence that is less reliable. We can use the `requests` library to download the file to the local directory from Python (or you could click on the links at the main URL above which is faster). The code below will download the file if it does not already exist in the directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "streaming-success",
   "metadata": {},
   "outputs": [],
   "source": [
    "# range-change file name locally\n",
    "a65_rng='ALOS2_A065_20180416-20190708_masked_range_change_m.grd'\n",
    "if os.path.exists(a65_rng) == False:   # if the file not here yet\n",
    "    # store direct download url of range-change dataset\n",
    "    url='https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/JL9YMS/PUWBZX'\n",
    "    if is_downloadable(url): # does URL point to a file?\n",
    "        # Download data into a requests variable with the above url\n",
    "        r = requests.get(url, allow_redirects=True)\n",
    "        # Write the content of above request variable to the file name\n",
    "        open(a65_rng, 'wb').write(r.content)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acknowledged-document",
   "metadata": {},
   "source": [
    "We also need the three components of the line-of-sight (LOS) vector: East, North, and Up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fabulous-attribute",
   "metadata": {},
   "outputs": [],
   "source": [
    "# East file name locally\n",
    "a65_east='ALOS2_A065_20180416-20190708_east_1asec.grd'\n",
    "if os.path.exists(a65_east) == False:   # if the file not here yet\n",
    "    # store direct download url of East dataset\n",
    "    url='https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/JL9YMS/CFBGNV'\n",
    "    # Download data into a requests variable with the above url\n",
    "    r = requests.get(url, allow_redirects=True)\n",
    "    # We are writing the content of above request to the file name\n",
    "    open(a65_east, 'wb').write(r.content)\n",
    "\n",
    "# North file name locally\n",
    "a65_north='ALOS2_A065_20180416-20190708_north_1asec.grd'\n",
    "if os.path.exists(a65_north) == False:   # if the file not here yet\n",
    "    # url of North dataset\n",
    "    url='https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/JL9YMS/PBIERE'\n",
    "    # Download data into a requests variable with the above url\n",
    "    r = requests.get(url, allow_redirects=True)\n",
    "    # Write the content of above request variable to the file name\n",
    "    open(a65_north, 'wb').write(r.content)\n",
    "\n",
    "# Up file name locally\n",
    "a65_up='ALOS2_A065_20180416-20190708_up_1asec.grd'\n",
    "if os.path.exists(a65_north) == False:   # if the file not here yet\n",
    "    # url of North dataset\n",
    "    url='https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/JL9YMS/YDEEMX'\n",
    "    # Download data into a requests variable with the above url\n",
    "    r = requests.get(url, allow_redirects=True)\n",
    "    # Write the content of above request variable to the file name\n",
    "    open(a65_up, 'wb').write(r.content)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "legislative-hospital",
   "metadata": {},
   "source": [
    "Now we are ready to create a CSI InSAR object with this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "resistant-dakota",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------\n",
      "---------------------------------\n",
      "Initialize InSAR data set A2_A065_coseismic\n",
      "reading from GMT grd files ALOS2_A065_20180416-20190708_masked_range_change_m.grd\n",
      "Read from file ALOS2_A065_20180416-20190708_masked_range_change_m.grd into data set A2_A065_coseismic\n"
     ]
    }
   ],
   "source": [
    "# Instantiate a insar object\n",
    "sar = insar('A2_A065_coseismic', utmzone=utmzone)\n",
    "print ('reading from GMT grd files',a65_rng) # print is a little redundant\n",
    "sar.read_from_grd( a65_rng, los=[a65_east, a65_north, a65_up])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "marine-shame",
   "metadata": {},
   "source": [
    "Check the dataset for NaN (no data in GMT) pixels and remove them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "collaborative-jacksonville",
   "metadata": {},
   "outputs": [],
   "source": [
    "# clean out NaNs if present\n",
    "sar.checkNaNs()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prospective-arkansas",
   "metadata": {},
   "source": [
    "We can use the `plot` method of the InSAR class to plot the dataset and check that it was read correctly. The `input_check` flag above controls whether this is done. This step can take a while for large images. The `norm` parameter controls the color scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "military-tobago",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/shapefile.py:391: UserWarning: Shapefile shape has invalid polygon: no exterior rings found (must have clockwise orientation); interpreting holes as exteriors.\n",
      "  warnings.warn('Shapefile shape has invalid polygon: no exterior rings found (must have clockwise orientation); interpreting holes as exteriors.')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAEHCAYAAAANnZUmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABeJElEQVR4nO2deZhcR3W331P33u6eTbtsS5Zt2dgYG4PN6oDZtwQnbAlJzM4TEgiLA/lYErYACQmBEAgECPsS1hBwgBiHJawxxhhsbGPjBS/yKkuWNNLs3fdWne+PUy21RjOjkTTTMyPX6+c+03errlsuna576tTviKqSSCQSie7gFroCiUQicU8iGd1EIpHoIsnoJhKJRBdJRjeRSCS6SDK6iUQi0UWS0U0kEokukoxuIpFIdJFkdBOHHSLyHBH5zjyUe7WIPGauy53iez4tIm+f7+9JLAzJ6C4CROSHIjIoIvVJx18rIleJyLCI3Cwirz2AMkVEbhKRX09xri4inxSRIRG5S0T+36TzmYi8XUTujN/9SxFZEc+9UES8iIx0bI85qAefJ1T186r6pHko976q+sO5LvdQiH3nTxe6HonZky90Be7piMhG4JHALuCpwH92ngaeD1wJ3Av4jojcpqpfmkXRjwKOAHIReYiq/rzj3FuBk4DjgKOAH4jIr1X1W/H824CHAw8DbgXuC0x03P9TVX3EgTxnIpEw0kh34Xk+cDHwaeAFnSdU9V2qepmqVqp6HfB14KxZlvuCeP0Fk8uN3/l3qjqoqtcAHwNeCCAiK4FXAX+mqreocZWqTnAIiMgxInKeiNwtIttF5APxuBORN4nILSKyVUT+XUSWx3MNEflcvH6niPxcRI6M55aLyCdEZLOI3BFH5lk890IRuTB+FhF5byx7l4hcKSKnxXOfFpEPicj/xBH7T0TkKBH5l/jmca2IPKDjGTaJyBPi50xE3iAiN8a3gUtF5JiDbJsHiMhlsZz/ABod51aKyPmx3Qbj5w3x3N9jP9gfiPVvt+n7ROS2+CZzqYg88mDqlZgfktFdeJ4PfD5uv902KpMREcH+gV29vwJFpBd4Zke554hILZ5bCawHrui45QpsNAtwP6ACnhldD9eLyMsnfcUDRGRbPPdmEZnxjSkaw/OBW4CNwNFAe7T+wrg9FjgB6Ac+EM+9AFgOHAOsBv4cGI/nPhPreSLwAOBJwFSv2U/CRv33BlYAfwxs7zj/R8CbgDVAE/gpcFnc/wrwnmke6/8BzwLOBpYBfwKMTd8KUxP/v3wN+CywCnvT+YOOSxzwKeyt5Fjs+T8AoKpvBP4PeIWq9qvqK+I9PwfOiOV9AfhPEWmQWByoatoWaAMeAZTAmrh/LfCX01z7Nsw41mdR7nOBuzH3UR3YCTwjnjsGUKDRcf0TgU3x87Pj+U8APcD9Y1lPjOdPAI7HjMH9gF8Dr99PfR7Wrs8U574HvKxj/+TYJjlmyC4C7j/pniMxA9nTcexZwA/i5xcCF8bPjwOuB34LcJPK+TTwsY79c4FrOvbvB+zs2N8EPCF+vg542hz0gUcBdwLScewi4O3TXH8GMNix/0PgT/fzHYPA6Qvd39NmWxrpLiwvAL6jqtvi/hfY1xWAiLwCGxH/rqo2Z1nul9XcEk3gvI5yR+LfZR3XLwOG4+f2SPJvVXVcVa/ERqVnA6jqTap6s6oGVf0V8LfYqHomjgFuUdVqinPrsRFwm1swg3skNvr7NvClOKn3LhEpsFFfAWyOboedwEcwH/ZeqOr3sZHhB4EtIvJREel89i0dn8en2O+f4ZlunObcbkTkwx0Tjm+Y4pL1wB0arWNkd3uISK+IfCS6X4aAHwMr2q6Uab7z1SJyTXSn7MTeFtbsr66J7pAm0hYIEenBXm0zEbkrHq5j/6BOV9Ur4nV/Avw18ChVvX0W5W7ARncPFZH2a2ov0BCRNaq6TUQ2A6cD343nT2eP2+LK+He2mp+KTfjNxG3AsSKST2F478SMaJtjMbfBlnjt24C3iU04XoCNMC/ARrprpjHke1dQ9f3A+0XkCODLwGuBN+/vvv1wGza5edV+vvvPMbfIdGwGjhYR6TC8x7LHoL8aG/2fqap3icgZwC/Z0+Z7/X+K/tu/Ah4PXK2qQUQG2f//o0SXSCPdhePpgAdOxV4ZzwBOwXx0zweLNwX+AXu1v2mW5T4Pe50+uaPcewO3Y6/gAP8OvClO0twH+DPsVRtVvTHW4Y1ioWWnYH7Q82OdntwxmXUfzHh9fT91ugQzLv8oIn1iE2TtCcEvAn8pIseLSH983v9Q1UpEHisi94ujuiHM7eBVdTPwHeCfRWSZ2GTcvUTk0ZO/WEQeIiJnxhHyKBaF4WfZljPxceDvROQkMe4vIqsPopyfYj8yfyEiuYj8PvDQjvMD2Ih7p4isAt4y6f4tmMun8/qK6M4Rkb9h77eaxAKTjO7C8QLgU6p6q6re1d6wV+HnxMmpt2MTSD/veEX98CzK/VBnmbHcD7PHxfAWbCR1C/Aj4J90T7gYmHE+Dptw+ibwZlX9Xjz3eOBKERnFRpznYYZyWlTVA0/BJr1uxX4A/jie/iTmRvgxcDNmFM+N547CJrOGgGtiXT8Xzz0fqGE+5cF43bopvn4ZFp0xGJ93O/Dumeo7S96DjZq/E+vX9oEfEKraAn4f80MPYu1yXscl/xLL3YZFuXxr7xJ4HzbpOSgi78fcMf+D/fDegrXnbQdar8T8IXu7khKJRCIxn6SRbiKRSHSRNJG2BImTJf8z1TlVnW62fV4RkWOxV/2pOFVVb+1mfRKJxUpyLyQSiUQXSe6FRCKR6CLJ6B5miKmHbRWRqzqO/aGYLGEQkQd3HC9E5DMi8qsYTP/6jnOPEZFfiMi74v7TRORrHedfLyI3dOw/RUS+0aV6P0dELu/YQoxfna7s14iIisiajmP/FJ/v0XH/v0Tk6R3nrxORN3XsfzWGcyUSh0QyuocfnwZ+Z9Kxq7CwpB9POv6H2LLi+wEPAl4SFyEAvBTTeshiPO5F2HLeNg8DhuKCAzBVsp90o95q0o1nqOoZWFzyJlW9fKpCxURonoiFqrWP3Sd+fBTQ1pW4KD4DMd52hH2f96IDfKZEYh+S0T3MUNUfAzsmHbtGTaVsn8uBvhgT3AO0sJhTsL6hQMB8/3cDu0TkxHj+aOCrREMV/x60UTrAenfyLGyBxXS8F3gde6/cyrDn6lxN9xP2fpbzgbVx4cPxwHiMd04kDolkdO/ZfAVbpbUZGwm+W1Xbhu/jmBF1avKPxP2Hi8jJwG+wYP2HR6N9f0zdqtv8MdMYXRF5KqZr0KmohqpejS2NvhD4t3j4UuA0MdWvh2Mrxa7DVgke6ig+kdhNChm7Z/NQbEnsemAl8H8i8r9R1Obb2OqmTtqjwQwzSpcAf4NJK16nh6i5e6CIyJnAmKruo38gJm/5RkzacR9U9dxJ+00RuRp4IKZI9i5see3DsedLroXEnJBGuvdsng18S1VLVd2KGdUHz3B92+/5cCx7xDAmuP0YFmYkeA7TuxbuhUlQXiEim4ANwGUictQM5V2E+XkHVHWQOJInjXQTc0gyuvdsbgUeF/2WfdgI79oZrv81Nip+JKZ0BXA5pqLV1ZGgiDhsInDK1EWq+itVPUJVN6rqRkzv4YH78cv+BHgJewTer8Ta5FhmIR6fSMyGZHQPM0Tki9ir/8kicruIvEhEniEit2Mz8N8Ukbbb4IOYXuxVmD/2U1E/d0qi9ODPgG2qWsbDP8Veww/J6B5gvcFGpLdPVl8TkY93hpcdIBdhz/JTgCgbuRX4haqGgywzkdiLtCItkUgkukga6SYSiUQXSUY3kUgkukgyunOMWMrwzWLpr68XkT/tONcrlvJ7m1j+qskrxKYq7yQRmRCRz+3v2kQisfhJcbpzzzuAF8W4z/sAPxSRX6rqpcBHsTY/BVt9dcYsyvsgC7PoIJFIzAPJ6M4xcbXT7t243UtERoCnAhtUtb3U9tKZyhKRc7D06RdhqW4SicQSJxndeUBEPoTlvOrB4lkvwIRbbsEy2z4PW3r7VlX96jRlLMPSmz8eeNEM3/Vi4MUAJ5544oPq9foh1V1VETm0xLGpjLkv4+qrr/62qk4WBEosQZLRnQdU9WUici4WX/oYLF34BuA0TCRmPXtiT3/doW3Qyd8Bn1DV22b6B6uqH8XcFpx22ml63nnnTXvtbBgcHGTlypWpjEVWxsknn7xm/1cllgJpIm2eUFWvqhdixvalWBrtEni7qrZU9UfAD5hCGyBqwz4BU8hKJBKHEWmkO//kmA7AgQh8PwbYCNwaR7n9mK7tqar6wLmuYCKR6B5ppDuHiMgRInKOiPSLSCYiv43pvX4fE+K+FXi9iOQichZmXCcreYG5C+6FRTecAXwY+Cbw2/P+EIlEYl5JRnduUcyVcDswCLwbeJWqfj1qFTwNOBvYBXwMeL6qXgsgIm8Qkf8BUNUxVb2rvWFZDCaikHgikVjCJPfCHBKN4qNnOH81e6eA6Tz3DzPc99ZDrlwikVgUpJFuIpFIdJFkdBOJRKKLJKM7DSLSEJFLROSKmAb8bR3nzo0puq+WmKJ80r3HiMgPxNKaXy0ir+w491YRuaMjffjZHecmpwXfGFOHn9txzQdE5IXz9uCJJUXqp0uP5NOdnibwOFUdEZECuDBOdPVgE2L3j/oKR0xxbwW8WlUvE5EB4FIR+a6q/jqef6+qvrvzBtk7LfingR/F/a3AK0XkI6ramtMnTBwOpH66xEgj3WlQYyTuFnFrRyf8o6o243Vbp7h3s6peFj8PA9dgKctnYqq04AB3A98DXnDwT5M4XEn9dOmRjO4MxFjby7Ff8e+q6s+AewOPFJGficiPROQh+yljI5ZN9mcdh18hIleKyCdFZCVMmxa8zT8CrxaRbC6eK3F4kfrp0iK5F2ZAVT1whoisAP5LRE7D2mwllrDwIcCXReQEnSLvkYj0Y1oLr+pQFvs3TFdB499/Bv4kft+5k8uIx28WkUuw7L0z1ZfBwcEDfs5OxsbGDun+VMb8lDETS62f3tNJRncWqOpOEfkh8DvYwofzYue9REQCsAZ7vUJMbPzxQB/2Kna+qp4Xz/UCbwH+CHsNvDZetxciUgc+hK1AWyciv4z7f4WtbJsSETlkYRUglbFIy9gfB9JP20Q/8FeBz7f7aSxrS8c1HwPOn2U1/gH4CjP003s6yb0wDSKyNo4cEJEeTIDmWuBrwOPi8XsDNWBbx63vwHQTvoZ1vkeLyIPiuY9iCmOnAKuwV7Srpvj6HLgNOAfzs70ZeA+wCfi9uXi+xOHBIfRTRESATwDXqOp7Jp1b17H7DKbup/sQV1j+mtRPpyWNdKdnHfCZ6J9ywJdV9XwRqQGfFJGrgBbwAlVVEVkPfFxVzxaRRwDPA64HVgNfFZG/x0TML8C0GBQzoi+Z/MWqOgq8NfrZiN97c7zv7MnXJ2bmpvXPpnfsVI7a+aaFrsp8cND9FDgL66e/ij5hgDeo6gXAu8TU7qbtpzPw95iOdGIKktGdBlW9EptYmHy8BTx3iuN3sscgPhuTcrw31vkexR4R883YaGQz8BlV3TxDHTYBp4nIkbGsC6YYkewWMV+/fv0BPePhxLW8G9n4bULWIBcFVbKbHsrNJwwx5E6FPhhc9eecctOHF7qqc8qh9NMoPTqlWLOqPu8A6rAJ04pu719BeouelmR054E5EjEHdvvcPo8Z6Gun+K69RMzn+FEWHVcd+RhafcfTK1vJK4fDU2Z1yHKCHIWgiIcMT3ANdrkcJ9Ys27KjFrj2iUQyuvNGnFG+UESey74i5hXwIxFpi5hPaXRFxAGfxV4PX9GVii8CbuJihlb/LzcffxmFQJCcHI8jp9CjudHdi5N8IMuboJBJgYpS4cipKGiieFwY4wi/ja2ZvQFsrG5a4CdLJJLR7QYHI2LeOclxJHB2lIY8rPjZyqfRs7yi5fqoZ56KAqUG5BS33I8sW4nE/zIqwNtIloCjQrRJFjyaKaJKjseT4V1JzTcpsxarqx9zdKWcIJ4KGClOIZTTvlgkEvNOMrpzSFxq+TgsvGYc890+C/PxdoqYvwM4E3M9vHaa4v4Ni3J4gqqOz2/N5w9F+eGG5xPqBf3iyaSkRqCkhz5diYjSSwApgAwFAjkWbQdKGZ2DI4AnD00qdyo1tpGFFgo0qlEynaAXx1FFwIuDDAaZ4EjxgJIBeRWoHQnbbl+QpkgkgGR055r28ssPYxMJtxBFzAFE5GnAx4G/juf2EjEHHqmqTxaR47DZ4iZwl+xJTPkSVf18F59n1lzOBxjaeBGaNei99f40N17FZk5gHZupy3JACAg9NBEUTy+OYUQqBKjwZEwwQYai1BkGdpLrOC1qQJOCXRSZ8gfhP2gC33XLebzfympVMoStkrG9coznnl51BBTyjKzZwgMI1MICNlIiQTK6c8pciZir6i1MM6u80NzkzmfTcV9kWdYEaeBxKHVa9OJkLZ4eRBqILCdQkONQBPDkVAhNciZo6AgZoxCEMer0yBgZTXIEpUXBsRyjmyFUoAEE6ijOwYhkXN0KnM1WekVoieBQVktgWIRbyTiFwLA4qAICZBqn05MUS2KBSUY3MWtudd/g2ntdzhHag3P9eGo4IKA0oiPAxrMBxwQBmIhj2F5GEMbxVAgVddlJowpkOsqGZpMel0NNCdRRlzOsgYwK7zyFV5wKkitjCtuCcGQG20NOT9kiD1DEnrwCz/panV2q5D6Qq8dVkAXbfA5wA3DigrVj4p5NiqVLzJodx32DURkgExu/OsYpGEUZxVPh2IVjJ54xHENs5Eq2Uqef2yjYSg+DHMUuNrKT42hxvBtlXauF0xyXgXrInKcuFU4CSkUWwIngBEIZuKjKcQFudBnrfAsCuACuBeM4tiA0W8KRQQnBk3kld0ABvmHPcfzK313Qdkzcs0kj3cSsKbK7OZbryBkE9XHSy9EjY0AB2mQ5TSo9lpIRcgLL2ckW4N44eiho0aLCU6cic0BvQV4GVEtyyRH1ZOIheArvyLIKB0gF9eB5qozziyxjPUIzy6mrR5wiAfqrQCODVt7imrzOCsArZB7wUI/PIYd9NHNiMZOMbmLWNKTJMm6hKRXjlGzAIzjGcCjQJwVj9OIEangKmvSS0SSjRoanRYOKGg4hJygU2jJx1gC586CChgwXFAklOQEndj6jwikcXcEvsoxhCs4KFbmaIf0VjtXBsRo4pSrZJpBl4EM0vJGwKL3liXsKyb2QmDVD5CjCcoSNgKfAU2O5CisUVFs0dYxtwFaUGkINZQPgYtiWRC+wUOHwiIDLMBcAgPdIcxx8RSO0cAqhrMhDifMB8YH1VDzFlxytnu1OqBQqhWUe1oWK3FdkqogIXoGqHYBmqICJtiUS3ScZ3cSsqalg0bQeR8adNNhKwYTkNEXIpKBXMtZRcipN6gTaQRhCabGyeHMXACEEQojRBZkjQym8JwMyVXypZK0WRSuQlYpGN4GUgZ0+MIqjrkIVwAP1PFDlNurNgicPSlbZ61zusN7uQDwc3/++BWjBRCIZ3XlDDiJh4GJP+NfUJi0cFSZodTIVGwnUUDJylJx+HBkWsKxxkYMHmvFviF1unIwhl9MMDg2ezFfUfEVZga9AAjQqJW9CVlnkWOahlcFQE4aqwECYYK0LDDi4Lof+DEQgBCgDNslWWSevgrkZ8Fa2HCYuhsOxnx3uJJ/uHDNJxHwL8E7gM5gOw/eBv8RUnsaBK4B3x/v2SvgnIr8CPhKPvUdEdqrqZ7v2IFMwMQa+v0UmFTsRGmTkBEotEGw0W1EQyPHUCHGUm7Hn9f56HEeSsRJPLoEeUWi10NI643Ix4zgiWCault1bAK4O2YRFKhQTsKwGF/bAMQKnewgKIzXoV6gr7GRPsLNgxluwkXBZ6167zTMHlJhycj8jJZbsOmmkO/e8A9ioqsuApwBvBx6K2Y2/Bh6ILQ1eBfxlR8LAyQn/PoiF8l+DafB+RETu28Xn2IeHb/4OIWS0FJarp6GeEhgTYVQcTSnwAp6AeQI844zQpLIVYQTW41mO5yoySrVlv8OZw4uNZLMKaiVkLRioYHkJtcrcA66Ko95gv1gt4LdacCzRc+BgoIxGPoaSQXvUDeSgOTRr8QZGWOocRGLKlFhygUlGd45R1avbHR3r1KuxznwpJhq9E7g/cBE2gn1I+z72JPz7FPAHWLaIALwKUyh7freeYzo0jFMPAfH2yt5SZYiC2ym4jRotoCKjFce+DXqokwGeQQIFQg14CBUrqMgksJqYs8jbKDaPixmKpjVggRnaoKBN6G3B+hb0NWGnh/Gmxfi2AxSqALc7CM4m6YKDkEOV2ZY5cB6OXvOgqR5xyXEgiSlTYsmFJ7kX5gER+RDwQuwV75fYiPe7wF2Ylu7rsGVR59GRMLCd8E9EHoDZkJthd8K/3wBPxIx153ftFjHfsGHDvCem1PwRjMVYW09BcDlryFilGR6lKRlh+/EE7L13DMUBvQQCUOIZjga4phWrq4rMVxTeRrcIOIVdw6dya2nDsWsyOFpgQKCSPT7frIBMYCSDIUfHtJ2xa+hUbovLf9vaC87bNQpIBjt3zdxeh2NiypRYcmFJRncemELEfCtwJ+brvRHLfVYC3wR2sG/CwH5g16Riv465Jy6e9F17iZjPdxLF69ZfwvJ6kzolGRmOgEoB5AxpBfTgNZAfdzENEZZTUhCifhi0cAwSWE2LFWEMLVusbLbobdriBWkCLbithGNWfoeyBVkPrAB6BEqBllgEQs1Bntvk2phAo4gxuO3wiDvhqA3fsYo7THu3tee0z0Hcvx5Se8yWxZqYcgpSYsl5JrkX5o9VWDK/DcBfAMcDFeYyeJSq/gi4BBhgUsJAzNm4bNKxMcwQL2jCv9+69Tt4D63Q9txCj05Q0xEKPKsYifG3Teq0KPBkVGR4SjxNAsuAATzDGo0fgEBZmQvBxeNOoN4DxwVYHiD3UFQw0bRfLFG7uVbBshK0hQ1hHXiBkNmKtBBXoGUOtGauBp9Dq+5Ycdxp+z7kEuJQElNORUosOf+kke78sQ6LWjgOe/M9H3tt+xTwsZgwcAPwuTga6eR67P/Nxo5jpwPfwtwWC0pLhhnwnkYVaAhscY6W5hzlR/EC/TRZxjhKhgsl9l6fU5BRANsJKBnDWQ/9ouyqWlzp9sw2irmDURdDwNr+AMxwHpGZj1fBWjaANMDl4G0FBllmv3DUbbVcBagqiqJFHGvkOW6fpl9yHFBiylmWmRJLziPJ6M4hk0TMr8bcAedhxvYC4BHAazDDeSZmRD8wuRxVHRWR8zBN3TNF5Cws/Ofhqvon8/8kM+OaFS4r6A8t6tri3pXZvjIHn9cY1xb1VgkSkFyokBijm6E4VgCOjBWMUXNCXm/wAD+BevilN3fCsgJ2FdCj0FMBEkfEAXAWqaC5xe9WmR3GgXdAZgY5F6hCIBfBI3gVqBdIZZHGVBWtvMDGzUX3G3IOONDElNOUsYmUWLJrpIadW9qhOrcDg1gM7qtU9esx3c7TsBjdXcDHmCRiHuMr27wMm4jbCnwReGmceV5wdHs/hR/DT7ToGYT6EBSj0LcDBkZa7GiV7KpG6Q27KMoRQmjSUmVcHcM0GKafEWqMU2ecglKEqrDR6f3q0B8TRzRyKAS25HB3ZoZWbBUx5Dahpjn8PIc7AlzfgmrC3BS34xjFwheaWUZVq1P19OKzgjIvaOZ1yryG00Bx3CMWuEUT9yTSSHcOmSsR87i/A3j6XNZvrjhu/CLynSfTAzQCjFewpYQTYsjXsSNw/DbPWAPGai1craJZr7HDgacXyOiJ/oKKGmSOUANxTdQHjq7gRgfjuYWPrfbm5x0Wi1aoO/ATmB84h4cE0Mz8tHcB24HjCWzJzM9QAUEyvC2xiMplLTKEQJ3KpbFHonsko5s4KLbvgo1ivtKBEoY8fFPhMS2oj1laHG0AAzDaG2hlo4wWaymlH6gxQUZFQYMWLaA/a6F4yCHLKooCbumFIoN1YiEfDypthkcDnBzMf9uMo+K6s4m1I+JI2ddrBBxDWY7icNTwZJQIjkAvDUvng0MkZ3jZU+gf+u+Fa9DEPYZkdBMHxbIWlBkMl7YYYX0Tji2BcdhVQZiwpU+9NaivB60roSgY1jpN6SPQSx8jLGOMFdRo4ejPWoxS0pM1aeUZ22s1slyp8pKTx6FVwn08UAdfWC7LepR9rAR8DZo1h89y1GUEcTjNKaWglIxAhtMKcIyJjXMLHC2pIWs9DC1smybuGSSjmzgo+jNwTetAAxPgRoFRYBiLHhiDbS3YfgQcVYdVPS2uqw2ztXECg6xlB8tYKXezjh2M0aSXXoYZYYAhJhCudT3cv6hRhcB6oEkJDv4Pi3IIHrTHvksqW/RQAq2g7CgKbpFejpCMpqtT4ggUtiBCAhNYMve76eMEPEqWNHYTXSMZ3cRBsbNpr/HFCBb9OYxNDw5jqxw2wZEBjrwF9HoI94L7nHELPz32vvzMPZLeImdFtpPNcjPr2MGRMsQKtrGSGqsY5RgKHL0sc03G62pi5HlJQ4UbvHKsgGaOAtN1DBmUDkayBr04TlRhkDpbGKAXwdPCk2Hr5oRhelhBQYllIm5JxsDCNWfiHkQyuomDorETagXIGOg4Zmx3AFtBB0CvAt1hWgm+MoHxB62H1edcwOue2csFW9/IwIYhVjRWcLL7DRP1HQxTp0WOp4d+enEM0EvBGC1cPsEKcpYhrMzVtB18RRmCOXQBT0ZNHE0cTRUqzalTo4VESXWPx9wMDQI5Ga14nw+Hj+xYYnGTjG7ioPDHXEe56WRqoxCGINsBbAG9DVgD41fBHVtNYOJqzPPATZBdCH9361f49dNPZ9Nlz6B2yijbj9hIoQ1KqUeFhhGUXrayhpUMs5oR+hHG8kAfSrMtFJk1yTMlsyk4LCYhp4nQop9RybiNPo5AbUSMKaC5qNLgyRgHlgNelmacbmLpkYxu4qApR0EnoNoO+d3AHcAWqOpw51Zb43zDpHs88N/vgtdf8WZe8o5nsH3nvfCNnMFGH+ONnJZkIFtwDDDIcgYIjOOpUexO8QOW9sfHzGt5TP7eJKOkxgQFTRo4aiwjY5wsuhEsdsGWHk8gOHoQWmRUErrZdIl7MClAcY4Rkc+JyGYRGRKR60XkT+Pxtjr/SMf25hnK2SgiF4jIoIjcFZX8F9WPZNgJMm6RCmEbVOOgI9AaM/X27TPcu/nb8MYHnga33MTOK09h07YTuGnsRG7UY7hBNzDMAJ4BtrOMIZazix6GaTBOH+P0MkwvYwwwTD+3sYKvs5pbWMYYvYzSQ4sGZczhNkIvO+mlRR9NaoxToPThECqghVJKnZ+seU13Gi5xj2ZR/SM+THgH8KKo1n8f4Ici8kv22KAVqlrNopwPYavR1mEiW9/FVqm9f+6rfHA0m1C2IG+ClFCMgy8tsmCCqH0wAznwd69/Mjd/5ul88pZ/4Pb8aNQFfCPjJG0yxJFY5rSMFr300yQ3Dyx7NLiVnQgrqXEDGStRVqN4aig54/RSUcOT06KJAnU8JRU5vbuzUrQAWc4sJGESiUMjGd05ZtJS3XbSgnsx88BvKo4HPqCqE8BdIvItYEEzR+zDLhOl0TJm2C1t1OuDGdxyP7cL4K+A53zxa9Re7jjPP5HNQ6cy3sq5P0PcpkfgRQlklIzQoocaLRPFgSgXKSjQwFEnZ5ycbTEJZo06Y/QyQU5A4gQaBFoIGRlQ4ihx9JGZSG8iMc8k98I0TJfwT0TeKiJ3iMjlcTt7inu/EPVLr8XmaT4OtIVqdoiIF5HtIvKrzvsnJQx8H/Ci6JJ4E/Bk4FuTEwaKyIvjPb84VAHzA6X2jOvQCVucoDWgiHKLkd5ZlrNzO6zMb2F8Vz/Vf96LbZtPZOvYaraFE7hRj+UO1rGV1dzFAHezgh2sYJgBhlnGVazgN6xmiDUMs4IJVjDECkZYRpMau+hlgn7GWcYgy7mbOsMMMEY/Y/QCdfrIsd+Mxtw30jxziP305I7zl0eX2Kv2d39KbHlopJHu9EyX8A/gvar67uluVNVni8jzMJ2FxwIvx0RrvgE8FZNsOQ0YUNULYMqEgW+J9wH8HfBZTCP18ZO+ay8R80N43oNivAd6emxUmxUWLkYGa7EIsjHM1TATV34cjvz4pRR39IC/DS64NxMPGOaGNcK9V4+wOQMvDfoYY4BR+mLqS41y5ErG9t35IByOgJDj6GOCVZjsrmcUR2CAwDaUnBpCg0AfkBGYWJpZag6ln14HnAGW8gebCv2vjkv2uT8ltjx0ktGdhqg9OlXCv9ne77F/AH8NjKnqNQAi8nux3FcAm0VkmaoOsW/CwG8D/4G9jfwc+8fxzkN/srllLMByNbEZ7YU8g76GPcgK7GEGMQM83cu7YInj3vqoZ/JPD30ft03cB3Qlu/p6ub3RQgbqZAwyIb2Ms4IGTWpUtAPM2t3Yxy3HpCMH6GWEPpp4mmSMAo4mY6yll+2siqFnTTwKBJZerO6h9tMOHg/cqKq37Oe6mRJb/gRLbPmxg/j+ewzJvTAD0yT8A3iFiFwpIp8Ukf3lYTkJmNyRXwF8P35eCfskDPw8cAzw71jnfhuWamWfV8SFRvqh2QdVH9AHsgbyfljbD6cAJ2PO6WPY/y/8jhvhEzteBZeMwXaBq9axeevR3Dq+gRvCWdweTmCz9nC7HskWjmILa9jOarbRxzb6GaaH8Ri9MEIf49QZZjnD9LCdGhPUGKOfcZRh1rCFXobIGKOHFr1MSMZnj3zfPLfY3DNH/fQc7G2sk33uT4ktD5000p2BaRL+/Rv2uq/x7z8T/bWTRMzHsVxV9wb+XkTOxAZ0H8bSq38A67xvad/fmTBQRG5mjwj1IDa6mFGSRVXnPTHlZMZPP4/qyg8z2Ae9R4AcA9uXn4p/GvgxaIzDmmBJ347EGmUCe5DJ0Q33BtacBU9aezGnLl9taX9vy2HXCqi16MlhWc8R9FGSSxl1ygKWLCIArfjZVpkV2wrqWjIcTmClbGGdK2lSx7E2TqspPTTpwSblBCUEt1cbLuHElNP208nELBNPBV7fcXja+1Niy0MjGd1Z0Jnwr9PHJSIfwwzs7ksxEfMPY28RO4Bfq+q/i8izsKR/R2A257vAH2LpexCRNwCPVNUnx7J+P5ZzCrbG4BdY6p9p305EpOtJFFeykiJ8h+Ul9O2E4k5rhBNu+g4M2VLgsSEYH4ddwXy87W0UcwfUMbX21WdD7Xvw3Hf9ks+9+ZN858RT4Y8yWLkd+j39/UeRFU3qMkaPjLBCdtErowRx9EigRo5Qohrw9HJcyLh+/Q7qEhiphJuyiv56k34msDwWgRzoYYwVKI6KwZDx5Grv51+KiSn3008n82TgMlXd0lHW7s+zuL+TlNhyPySjOw0ishYoY0duJ/x7p4isU9XN8bJnYMkngX1FzEXkS5hvFlX9IvDFzvtF5C/b908hYn65iJwDnK+qp8Xrv4wlDLxkPp75YPEKVQ4jNVhZx7I6LMeEyWs20VYbh4ExGBmx/GYtNcM7jo14a0C4APyZ8NDP3M3q+/8rP3r5EM3r7gVX3AdqI4ysvBqq4y2lxPJeho5wLO93LC8UJaeJqYcFKkZCYMV4YHu5kWbWIM/G0JqjUs+g9pMzwQpKgmvRR06GEoB+WVoet4Ppp1PwLCa5Fg7w/t2o6rUi0k5suaj66WIhGd3pmS7h32dF5AxsQLcJy2OGiKwHPq6qZ8f9XuCJ7fMdvGuq+2fJokwYOJZDUcJyB1WvhY+5lSADIHWbYHOjJr+4ogKGLNvEwC7YPGFG9yZgOIOjfgdWroIjJuArX38Ol/3ZA3nL1/8a8rNh7Akw2LKpnHKYwV6l3m+SNhMKGvOcZYzSapa4fJSydSyjE3ewJs8IawoqJlglI2yVguUoGUJFnSaenIJBavyGmziJExawRQ+I1E+XGMnoTsMMCf+eN831d9Ix0aWqY8Dq2d4/TZmbWAIJA+88/Q2suPIfCMttQddQH5QrbJVa5sHVQZdhy75GgeXQOwplDdxt5ug+CZCjof9ixx1PDYRBOH4Qbs3hvX/0z/zgrs/xjas/Ats87FwFtzcg72F0udDykNcA9VEOvYfBrcLE8ATNdVtZv3qCtQyzpeplND8C0V6Ws5OaWDzDGP20aLGSFi2Bz228nLdtWhpGN/XTpUdqmMQhc8LqFzDYC1UDWA49DahWQrYS3HpgLcgKkGXAUcAq+9s8FY6/HzzqJOh/IAy/CEIeOOnTsPZ2GPgIPKF+GRf+4qEMHfVkTnvA/8CKzdDvYGUTelYyPLyS4bFexpoFzdBgPNSoXMayDZ6BWouhoYLC30Ut30nhdlJ6JZTD+LKO1waZQg+7UJQxehEyXL7k07InFjFppJuYG1ZBFqCvgOE69CwHmnGlWh2kD5s1a9myYR8gGwc9CvJtsKaA7ErLe1asB06Eda8CWQUvuuyDvOmk91NfNcb6Z97FndeNwNjx5kD2SjmWkxeQZU2cqyiAWi7Q38vA8pJbxk9hK0KjHqiHJhOuTlMyWjrGEVJQo06TghY5q1BSxFNiPkkj3cTckMHgCnMt5DmU/VANgPTGEe4yCDGeV5dB6IViNWRHAveBxgnmIhg9Cfy9gdUgR4LrgUf/b+Adq1/BuZ9/NX888W7oWQVrtkF+G0xYSuAQAD9OJp5cbGoulxp9/TUGlhe0anV2yQDDeYOqVaen2aLmClSUiZgnuIFHqdHHGB9Y9eUFbMzE4Uwyuok5oTr+QooeqK20iIUcKOKyYO23v20Xw0QfDK4GlwMN8A4Gc+hdBSc3bUIu85ZnjRuhkcEjngq//9Mm53zvJzyy74usW3kHhB6oMtBAWSo+5IBSqceFEiVQtYRd4wPU1VMgQC8TPavxvQME+tjCAJ46EyyjxTh1higkwGq4nCsWskkThynJ6CbmhAFZS7MG44VtrAQaZoBDbvYx1CHrh3wAVjUsrGxsFTSPgmIF3CZRim3YUgC1bgcdAu2D+n2gXAd9v6j42+a/cET9O5Cp9eBmP6GV44MjqOJxlFIwgUNyYXlviyrUaVU1Ms3JtEnQnKZm1LSPCV0V17PmeAocNZx4LjsxGd3E3JOM7hwjcyRiHu85R0SuEZFREblRRB7Znac4ODxQr0FWg1Y/VMthotf+Zj0w1A8jDQsho8/8vPUCGgFqw7B+B6y6E8KtwJ2QXwPVZVB+H6ovwPDn4LY/6OWo4Ql+8OMP8s1LzmDtym9AeTsM1hgf72d4vIdmVaPlc0pfMOyXsW1kGdtHVzFeNagECmkiklOxmlHpZ5tkDFHRRxEVdj2C31syLZGYI9JE2twzJyLmIvJETODmj7Eg83XzVeG5YlsvLNsJrSaMN6AnwFAOuQACy6POrgBUFlImmU3AZRNQ3wV+EMa2wEgLLrvYVpZ4LNJMLnoc97nvDay7+Va4Hs78QuD7l7yY7DQ440v/Rev6hzC+vKBv4xi9/S0G8hY9RZNlhVJpP4qSkbFDjqbJOCtkgn5K6rTox9NLgJiBLSOAKBet+T9+1z91oZo0cRiSjO4cM4ci5m8D/lZVL477d8xB9eaVI1dch2+dDGKGNs9hbYxYAMsK7CpwAchM+Nx7yzQR6oCHbDuM3ArbToLr7BDnv/9PWHbuKWzUS9DRaywm936w8umw6veBXTD6tWfwvdPhdzZcapoOzkTKV0QL72QnLRzKMD0I/bToI9AbFcsa+JhfwpQbcrx9crESicQckdwL84CIfEhExjAR883ABR2nbxGR20XkUyKyZpr7M+DBwFoRuSFe/4G4zHPytQsmYj4VowU4gbEajObYz3oD6APXB241+BXgC/ANczG4PrjhKBgKQB3GvYnEPgy4+OsvZe25R3NfLuf+cgtH9ZSEAdBjwD8d9PdAnwGbnwhDD1rGW459LT0DtzKsHmjSoIUDWtTw9KIMULfqxNQ/IeYF9gglGZ6MkoyAo8Ilg5uYY5LRnQdU9WXAAPBI4DzMhmwDHoKJ1jwonv/8NEUciemiPjOWcQa26uhNU3zXR1X1war64G6IquyPqse0dXfltjy46oFBB0MCzRi5MFzA2HLQ1aBrIayBe9dBT4HRB8LGZ8KGB8GXvns/Tn1qxW/xK07mRo5hO8szJV9eY/wI2Hk0bD0atm2A5avgt4aHeMjdl/La8v2s1itATdImRu7GT546HqGioMRGtaD43QpleVTXdRDTtScSc0dyL8wTHSLmzwVeqqrvx5TCALaIyGQR807G499/7RDHeQ9mdN/YheofNHc3QAq4d59NqlHCiob5bqlAm1A4GAngJ2wibXAt9Jaw2kGzH25fC9fU+jnxvuuA2ykYp4dxGqoEGrSKGttrHgcEFDcwwrKxcQZKeIDbxbKJn7Cm/25+zOspwyo8woQ4MhEGIMrbgFBS0KROoE5BRTNqMWRRAMcno5uYc5LRnX9yzKc7mfbUuOxzQnVQRG7vuGbJsIHr2FmcS1FA5iyKoVECIfptG+bTLZu2WCJzMOCg5WBwALaugc0Da8g2L6efrTgyarQQWjgpyMgZx9zENhvpyLNexgcq7kK5i4ItHElO4Gi9jo/ffAJ9tVU8aP0W+jKNCSlbZEzQwwQNmtQBF1NdQgtHExf1dXc7pBOJOSIZ3QNERBqYVmgUMeQrqvoWEXkrpsTUxPRy34ClDnsW8OwOEfNvY6qG64BhVd0Vy10PfA4YBp6D6ey+Q0Q+DJwOvAo4X0RGVLW/Kw97kIz3muugbEKrFcVoPCCQCeBgVQt2lTCSQdEPYw0Y7XNM1HtpuDr9tFCGKHFARoGjYI/weUlGizoeR04PYyi7ELaxnF+HdWzhWE6So3jCxu3U3BWcwAR91BmmB1PvtZk7weNo0aDCRdeD2z1zFnDOT59naBGwn/74Z1gaHYA3tPPxTbp/E9bnPFCp6oPj8b36Y8zB9lbgdcBGVd0ar1v0/XGxkYzugTNTIsAPY5kjTge+gKXpeZWqfr1DxPxYLK3Kf2MduM0XMR/7x7GMEX+HiUuvw2TyPolJ5v3l/D7eoTPcyBgZgJ4CpAXNcSga7IkGiKnIapg9a9ZgvGEGN6e+e6JBKSnIKahwCBkWe2a50OoEzGyO00dJwRZybmI1v3GnsoONrJXt3Nttw9OLp0AZ5giaeMbZk+qrpEARJP5jUFxMHu8JyOJPy37QiSk7eKyqbpt07C+Ac4ETsP744Xh8G/Bq4K8Over3TJLRPUD2kwhwVFUfPc19bRHzTcCDp+jkl2CG9QRAVLUUkW8C3wReCLxVVSdE9vFGLDrWbj+Xlcd8G+lVvCrqoIypDDWLbgbM5E3kUNVqVFlj94t/Fn2pNUCp0JjVATIEF50AnjEylIwx6kxQY5hemqxkDRV97GIZBYM0OYqSjCa9tCgAYRc1AgI0kTi+3ZNvMaNCycmoyBZ5vNgcJqacTLtBAnu7wD4JvFBE3qmqO+bge+5xpOiFg+AQEwEq8B0RuVREXtxx/APAR4A/x17r2oxgHf2Vc/oQ88oqmnlOM8+Z6HeM1W1l2uhyGF3uGOspGK85mnVHq15jyNUYkwYtCloKleagkNEip4XgUZQGFQUtoMRTkiEIjhyN+c5gLeMcw05O5jrWsY0T2MpR3Mkx3MXR3Mqx3MSRDLKMEQrGyamoM0Yvu+hhmAbD9LCLXnZSZ5xCxjn/mNcsdIPOSOqPS4s00j0IDjER4FmqemdMYvldEblWVX8cU18/apqvfD9wuYj8837q1fXElNOVcfctjycgSPBMVEqRic2a7f6dDwQVWmQ0XQ1BERxBxCJkd9yLIBqvtqCuEUDRmJIno6CGQ3AUZDga1DiCHA8UKMX2HjzrKagoqGIomMaoB0cDF825ReUGBLF4CCzVpSPffiRO3SG16yJPTDlv/TExNcnoHgIHkwgwKvejqltF5L+Ah7KfJH7xe74AvGym6xYiMeV0NMJWlh3/KwCqsqSkwNdqZLrnRb6JYyxr0FC1iAZXp0UNJ6AIjY0XITh8NILt2AJHRoEwTg0XXQ5KgcRpsBotCjwFp5NvvJQaJTkVe1RylYpAE0dJQUUDducHzoGMMtYvABPHXcWFx1zLU25670G3x2JNTDmf/TExNcm9cICIyNo4oqAjEeC1ItKpjTBlIj8R6RORgfZn4ElTXTcN78GiI5bGD2X4At6ZR3RMG2yXGpTgQ4UPlTkJRFDNQHNKyVA1H6poHj23RTS4tk7M4mczMiQGctnEWo2KAcZZzTbWs5W17GQ5LWqU9DNCL2P00qKHceo06WEXDYZoMEKdcXoZoZfReH6COqM0aJLHpRW5jNOXTSxcW85A6o9Lj9RgB86hJAI8Env9A2v7L6jqt2bzpaq6LY5EFn30QpsAqMup1W1Bwph4+kNOBoy5jCB1cvW28FaEUhwOqAgWvwTRfRDigt0QTa+jh4IWFTXG4+oyII6IbSw9QjMGlZlLwfQUKmz5b50KRyAnj+4Ej8czQUXA4ZCYrjLQRFmzeON1U39cYohNfiYOB0477TQ977zzDqmMwcHBQ34VbpdRnvgwel0LD+wqM+qU9OUZlQitKClTUuCrQEkOeQFkKI7xTWdS23gxgoWWZVQEZLe/VaOTIqOMBhdKcgJZ9M4qo5vOpG/jxShQUO72BxO9uDZ6rhFwMT6i/beIhjujtekMqo2X09RAfUg5bcsXFqRNTz755EvbMbSJpU0a6SbmjVyFZhxrNjp7muZ4UVrkFiGb59TikgSbzDKnQgOoaBJ1cMgwP6tFRXlyyjjdBbeSUaNiJRU+LuPdEwJWxvUNbdE3obLcFkhUWwiAowZxtVqFcB099BNYjWdclOZAL2zpRsslDmeST3eOmUsR83jfSSIyISKf29+1iw3vsygikyECWcjwmuPVUSOnhpCTkVPszlHWQ0nPbpmail4y+sl2T4LZWLaFo6IixIBUzzF41u4en1bYAqsM0Ogu8AQC5uiw77Tg07ZxFjTGB5tMTsl9GWMlnhYlFYFeGSeROFTSSHfumRMR8w4+CPx8rivZDXR7xuajoCaOGo5GVsalBiHGErTdBEYW15q1HV7tv548BncFAkK127drrgePw6FxxZqP5lOZoCKPy3v3RPe3jbGZ2nFkt8NC8ORU8YwJQjbwlDhWE9i59KQwEouQNNKdY1T1alVttnfZI2J+wIjIOZhew/fmpnbdpRj+IavxDGigph6njj1TCH630ldOFbVtzXGQQ/Su+hjA1UTQqL8gNOLI2Ba2mV/XwsJsMUWdkiwa0YwyRu56FE+F4uLotiDQgyJxQXFGxQjKDhw1PLVoxF2cUGssZhGGxJIhGd154FBFzGMZy4C/xda5z/Rdi0rEfDJChdAkr8ZxfsxCxmLYVx7dBG3xcGhry1isQi0aS2irgBktzGXRdgrUaUX5myp6kH10MWg03ebbdZh7IuxeCmFbQUUtrmrro+QommT4WD+NkQ6eJiU/WvtHXWy9xOFIMrrzwByImIOtIvqEqt62n+9aVCLmk5EQqNQRMo+4QO5KHE1CWUFZUZSe2m6NgxYtbLKsisv+286H0GFMe2I2CNvaIWEVQ3gymuQ0gVYUa/QUeHICRUzJU6OKDogKKKOETglxabE5PiqgQikJMdaiITm1ZRmJxKGQjO48oapeVS8ENmAi5iOq+gtVrVR1C/AK4ElxRLsXMb7yCcDBL4FaJMi4o5eSgoKWFDixaYTgMkKWQ9aM/lpbhWbxsm2JG4Xodmh7fDX6gdsv+m2TK2Qsi1eXcSrMEaI/N+BokdMkiwbZUSIoNYiqYiVRBw1zRdiqtAyNPwoQaNJwizZeN7FESBNp888Bi5gDjwE2ArfGwPV+IBORU1X1gfNQx3nD3XkJ4aSHIpJRw+NxCIFaFijVEeil0vYCXIuvLQQmdi/LhZwqGtccHxc5xNIZJ+fXCBspWRPVds05kUU1h7bB3DO+sDVu5n4o8QgZE+RRdUGi0GMzxjeY7Lmjopc64yzOlWmJpUMyunNIFA15HLbOfRwbrU4WMf8NsBITDflhW8R8Eh8FvtSx/xrMCL90vuo+nxQKXoC40gtyE50RKNXHaS/LW9aPUGrGhFrnrEUDm8Uxqsk7+ri4AUZxjCExZkFitgcf4281LnhQinhMoyE1zJTbVUDU8jVpR4nfRlR/EISKOsm9kDg0knthblHMMN4ODALvJoqYYzq538KU+K/C/LzPat8oIm9oi0+r6piq3tXeMDm9CVW9myVIeWeOU4t/LWJYVhXVvXaJsDlk1AL0hMC4Bw2OAQmspKRGkzpNcsYpaO0eJVQx4mAA5QyEZTha1AjxCgswM2lyQdFQgXraSXjaKO1EPdlu6ZzOUbEluhBq5BaSJiWXHfOMLrVc4nAkjXTnkGgUZxQxn+Hef5jh3FsPuXILydhFiD7EMlZSxlGo/UL141gewEkLqFPLBMTFOFvTDstRaihKkwmgwqQgc0rqTERTmeEJ+Cjh2AOM4+MkWYa4GiGOpomjYDAJSMtL4WP8QhlDyiR6l6vozhhnAodqQb023P02TBw2pJFuoitM3PCfZKGJD3C9OryaBHkeIMtK1GWUmTch8xhDUJFxBxnDMQKXuEgii5EJKOQaqGsFWoGWqCqVZlTq40SZ4miZV1gDLiiEFvgWVSipguJUEA3RIFtKykw9LgBBaZYV6idoeE+ftuiRA1nbkkjsTTK60yAiDRG5RESuEJGrReRt8fhbReQOEbk8bmdPce8xIvIDEbkm3vvKjnPT3i8i/xRjbh8d99tLh8/tuOYDIvLCeX34eSBjI2NSp+karJfoYQ0VdWnhJOBdBiqo7FFJaFBjPbWopGBugXZIV66BLMBoCb7yFN6TBUfhofCeGgqqoIGgNUAQcYgGMg049WTR2NaiqSdG7irBpCadA1fQKCBzjpp4eqSFLqKMSamfLj2Se2F6DiXhXwW8WlUvi3qll4rId1X119PdH5cMg6n1fxr4UdzfCrxSRD6iqks6XqkkUCmICg5nYWPtgDBVRBxO9+iIhQ79hLYyWEYgqKLqCOrRqmRHM7DcNWk0crwTVBwhZNGDoOTi8epM10yhVMjEkWuJ4ghBCa4gw8cYChen9tr+XqGQQEsqmlgdN/EDNvLYhWrKTlI/XWKkke40qHFQCf9UdbOqXhY/DwPXAEfv57Y9mRH3DiO7G1sG/ILZ135x4jRQB5xYFl4fJnavOQsI2rIFExoFazJawAQZIfpcTb4moyKXwHCZcdeOJmNDJWSBzFdkrQrXapL7CicWp+tCSR5aaPAElJoDRFFVgmCxCyoQfMwVYXG/OS4uwGirPNSBgkqEwY2LI1NN6qdLj2R0Z0AOLeFfu4yNwAOAn3Uc3ud+Vb0a6AUuxPJbdfKPwKujUPWSJdcmpSqot8kqVwMF5z0SmpAF1HlUW3GZbkWBksVluZkvKYKnphU1P8JA2In09bJmZUG9UccXOVLLcUWOZqbgQBBQ2R2iFtRRBshFIDP1M8RRU09dm4g2cepRlajQYFkrbLM0mIEaWbZ4BnOpny4tknthBg4x4R8i0g98FQsbG4qHp71fVc+dqhxVvVlELgGevZ/6LprElFPhb3kiTkCUqBgGPi5BcK6CIFGXIdDccVIUbpSol+BA215XcAqDLWgQGKvnjG+3ZcOi0B6Ajew4mTJ+LgiItiXMLXZ3T2lYQkxRvOZRr8zieSd2nMBu+UclxkoINZVZt/UiT0zZ9X56TycZ3VlwMAn/on/tq8DnVXV3Ooe4BHi/90/BPwBfYYakgbKIElNOVUa19rtkvVlcYAtoZol2xITLVcFNjKP1HnaKMLDR3IUOyLVCgyCUqDg0ZCxTbwPZvIZ4IC7tzVCCCIiwfOMPY5RERR4qVBVRoRLBidBCEAEVoRKHio2KS7FoCZHAso0/BSpTGgs5gxIYo8bK6m2H1B5zzVLpp/d0knthGuQgE/5JFDEHRrHMqkPxeHuGd1SiiDnwocn3x2vrIvIJ7BXuFDE93uOBXwO/N8eP2jXyzVcQQsAUDTKqmNVLCRYR4KDsqaPO7dYF07hGzIvgM1ApEM1wAj4rUFcDQDMgc2QqBBUKtUiFmqrlWNOAdyDBVqDVNSDBk4UAAsE5MpeRicNLgYuOhCz6kC0S2KEuUIilfr9s2VMWrC3bHGw/jdcL8AngGlV9z6Rz+71/KlT1WpZ4P51vktGdnnXAD0TkSkxE/Luqej7wLhH5VTz+WGJiPhFZLyIXYCLmz8YmNLYC/yYi12F6CgD/BdwM3ITNHk+V2C8HbgPOwSY33gx8GfgUJqCzZMmzHFwBzl7fFRuh26t+VNhVmzTLNIAGvHqqtkKuiwbWQSEBkYBTR66OnAzE4cSMdhAXR6/gnJnNqsij3KNNmoXc4cUhKmjI8eS7l26omrvCaWWhZsHjQkUWWvT4Jn1rRxeqGTs52H4KcBbwPOBxU4SGTXn/LPl7lng/nU+Se2EaVPVKbGJh8vHnTXP9nUC7w16NDSROBn6IGc1L4rkX7i9zhKqOAm+Nu6dhhd0M9Krqkv6hDFoRNCNoQCWHUEZjCc4BXs1IBouzdRKioliGi6FcxNAzh5I7bNEDgFaIs/wPGiUWvAg+QD3qLji1TBDgyLMMnCMTy8rmaRGdGbaSTRxePCG0RW6UWlyvVpKRu4UXvzmUfhpV8KaMOp7u/mmu3UTsp3H/CtKAblpSw8wDMgci5pPKOxK4N2bMJ59b1CLmkwlqmguV2ASa5gU+FzQu/yVGEyCKiqKS4xDyYCNPjbI2dZdRc5bvLBOQ0KTyJT4oBLVJN83Jyam3R9Yu4BDqGdQzAWd6C5m2oxeUunqK0ESCrXBzKggOVAlq5tg5Ty0P9KRMEomDIBndeWCORMyB3RMdnwc+E/1lk79rUYuY78NNT6IScBVoGcgrT6gqCIqEQOacTWzhIFiO4ELMIeBU8Ooo1eG1QrUCDYwpDGZ1MsnIQknZDk0jxFQUDqiTSQ2X1yBvUOUFwUkcEXtaYosuUFMVqxFsdKwVFR7NM4o8o8gEJ2r5iCVw5YanLWRrJpYgyejOE4ciYt5GRBzwWSxDzSu6UvF5JuOfCChFXpFllqHXCfGVPacURxVgTIVARfCBKggtpwSx+N5MW3jNaYZAJUpdAisAl9XQooe+vE7RTjApFUzyyJjIY4ZGnQcPZMEWQhA8WRAyEaAglxoDrpc+rVOnToajHiXWc3KyesoQnDgwkk93/jkYEfPOmeUjgbNVtZyf6nWfTHO85JQOMsbJVFH1SPA4MlMJ0wq8IkWMqQ2W1LJCKcz5i7ocm/CSGPVbi7G/IGKRuJZ9ogWa48hpL2kIksWk7ELpAXU4KcmzDA0thIJMKpCAlyYmZl5jzz8Zk0FvpLTsiQMkGd05ROZOxBwsOP0U4Amqepj9yy7JRHE+gCvw6vHiyDNPi4BToXJCyNumLTORG/GgigaLZJAMvMYElaqItkBycifg4yKMEAhZO7l7W78sj+kqBajRzAocTWrUAWU0M6nyesdqWkvkUwEFRcfC4DItvkocIMm9MLfMiYi5iBwHvAQ4A7irHdcrIs/p4rPMH+O28KDKa7TEk2dCEcPGnCrqbNVXoaChwHuHBlO2rTyUIaAioI5SwEuODyYXGYLQBJqZi7G3IN6TBQ/qCCpU6gnazrZmqdZz6mhMSdkPNGKiTFEPAXz096KOCkueqRS7RdMTidmSeswcMlci5qp6C9O4HQ4H6ndezMSJZ4F4Msl3z3W5ULeJKy9IUBtZSg3iaLISxeU1VD2qHqSOBI+ox4sFKGRSkStUFLvFamriUQ0mXy4CYtl+20uAs2hsIeyOR8ii39eJMwkHbY92S2y0a1f3EfjZwB9y5vB/drcRE0uWZHQTC4ITRxBicnRBsRhb52pkKFWWMyz1mIk3kDttZ1bDi8OrQ9TcBeIcmWIREAg4yCltlCoOyyFs6dXBUUX3QAa0YjJLM7wFAZus87vHvXYdEifcFLx6eoASxy4V+tdO2PtLIjELktFNLAgSAmSKKFTEZJUSYkZey1RWKwM1cbaALdio1JOBC9TU46Mur6Kmy6Ae56J4udgf5ycIkkEorHwNZKKoy1E1scYKsSXHYiPjijKuS3NxNJyBLeGgpQGhGY826JOc0i3++OjE4iEZ3cSCUI0poV+iwTXjRzS2qKeGI69ZlIHJ3MRwMQJohg8BR0AzMbFzZxklgkAIkKvF/fpKkVwhq6gkB2yxQzDJspjc3RZtEAQvjkxsqiNE7QdtJ7tUwYnJyQax2vbSYihlCE4cAGkiLbEg9Gy+yLymqjhVSs1QNb3bTIKphIW6rW8I5hQoUZqYYIVmNTSr20hUcnAZwWWoOLJMyDQgTsiK6KsNDlcFsqA49RBH1xU5PmRAiadlkRShotKwewTsFXzwCOMoTRTI0Ziw0tIMJRKzJY10EwtGQU5lKxJoOxkcQkUN1CQXLXVOAK1Tc5ZCnfjSXyGIOgoFxO1Ovi4BgivsKnVRKcybX1aEKhpSQRDJQUw2EhRPk0wzOwaEUGFj6YBqhiNHqfBq/l3VwPE3XbQwDZhYkiSjm1gwyijRCOCJCx5ULLMvSiVQqEUJ5OKi0KLgvEMQagQ0s9hZgqNSj1MImZJ7jxlRjxCoYa4E1QyvNXvFc20lB8gcVFKQUXF3K2PYB9aEEgmeqiwpyAjbQEZ+urv+OW23SCIxe5J7YZ6Qg8jSKvewLKvjgZitNxDUoZpBjGJQiUtxfcCVnlC1LIeZzY6hrqRyHlFLQCkSQFo416IIavo2qAnUtNPttErC2ATl6Dg6oQS1tOxKjaA5uVfcBBy16dOcfNvPWH3HZcj2N+E2XYpuugQZuWR/j9R1Uj9beqSR7hwjIp8DHg/0AVuAdwKfAS4Rkb8BSizMfjnwTlW9IN63V5ZVEfkV8JF47D0islNVP9u9J5l/GresYvCEHewIjgYZy/H0uJ6OdWAVIVeGtWBlmMBpFdeRKW63Bq/5e4PYqFl9ZvKQru2ssIUWQQN50SA42NWqaJUVx93xy4V69LnkgLIBT+5npGy+XScZ3bnnHcCLVLUZO/gPMSHyIp7/R2BoitTYk7OsfhATurkGWzr8ERG5bP6r3z36/H+CPol+B+BxYQLEmqmdNUIEBiSgrkFQb/G4AohEn2qJSoZT021wrsJh+riq0LpjACb2zhyzbnJFljCqqsCBZAOeKZvvT7Bsvh+b+5om2iT3whyjqlerajPuCrAaS039fx2X7S/L6qeAPwDeg/0DeRU2Qn5+Vx6ii3hpMiEjiLZoSU4VNIZxBcqgbPUZO0qHVhVOSzIxTdtMA+oDGhzSVOQ3H0N/cwn+N5dR3ngp4cZLCdvetI/BPRyRA8gGnLL5LjxppDsPiMiHgBcCPcAvgadgyfoA/gwTwroUE8DZJ8uqiDwgXnNzPH6ziPwGeCKw11S5iLwYeDHAhg0bFnU24KkYv/ksclGGkZgq3SIUmjtOwoVANj6OOGFbvUFeBoqdr5+htL2fvdvPMp9lzMSBZgNO2XwXlmR05wFVfVmclHgYlhttK/B94GfAa7DR7wcxEZzVUxTRD0xWH/s68NfAxZO+66PARwFOO+00XczZgKe8rvoHdp3wCILk1FAmQoUO1qnf/TpWViv3OGUq7L3hAKu21NrjUDiYbMBTkLL5zjPJ6M4fqzA1secCf4FNXLwTWKuqm0XkFVgqn69Oce8IMFncfAwzxL/HnnxrhwXLb7pw9+e++HeQQxux31MQkbVAGQ1uOxvwO0VknapujpcdUDZfEWln8z2s+tliIRnd+WMdFrVwHOaXfb+qni8inxWRM2D32tE3THHv9dj/m40dx07HpCFfOE/1TSxN1gGfiX5YB3x5Uj9TYBMmFTpb/h5ziyXmgWR055BJIuZXY+6A8zAf2V0xO/ALsJfkDwFbVPX6yeWo6qiInIf9QzlTRM4CngY8XFX/pCsPk1gSHGg24GnK2ETK5ts1UsPOLXMiYh55GTYRtxXT4X1pnHlOJBJLmDTSnUPmSsQ87u8Anj6X9UskEgtPGukmEolEF0lGN5FIJLpIMrqJRCLRRZLRTSQSiS6SjG4ikUh0kWR0E4lEooskozvHiMjnRGSziAyJyPUi8qdTXPOWKBr9hBnK2SgiF4jIoIjcFUWlU4hfIrHESUZ37nkHsFFVlwFPBd4uIg9qnxSRewHPxHQXZuJD2MKIdcAZWPzvy+ajwolEonskozvHTNLT1bjdq+OSDwB/hQmUz8Tx2Dr6CVW9C1vNdt+5rm8ikeguyejOAyLyIREZA67FRrTtlDx/CLTaKXr2w/uAc0SkV0SOBp6MGd5EIrGESUZ3GqZL+Ndx/jXRL7tmitvfhymF3YC5B+4UkddhWqU3i8gdwHrgwzMkDPwRJmQyimk5/AL42uSEgSLy4njPLw5VwDyx9DiUfioiJ3ckrrw8zkO8Kp5LiS3niWR0p6ed8O90zKf6OyLyWwAicgyWxeHWqW5U1etU9QxVPQn4Rjx8IvBZLFvEe4E7gT+fJjHly4FvYyPbrVgGiVWYHu/k7/qoqj5YVR/cDaHsxKJjLvrpGcCDMM3m/+q45L3t8zP00zbtxJa1OXuyw5RkdKdBjekS/r0XeF3H/kwci4mSPxQTM38N8DfAMcCXReSv4nWdCQPr8fy/YwkDvwvcBpxNItHBHPbTxwM3quot+7lupsSW38OkSxMzkIzuDEyV8E9EngrcETVHJ19/hIicIyL98d7fxpT8/xvr1KcBH8YyQHjgcuDzsE/CwPdho9vnxqI/iCWlvHK+njWxdDnQfjoN57CvCl5KbDkPpLjPGZgi4d/9gTcCT5ruFkxP98PYD9qtwATwN6q6HcwfhkUv3IhNsv0t8Cci8gbgRFV9ULzu92M5p2D51bZjKbJPmaG+Sy4xZSrj0DmIfroX0SXwVKAz62dKbDlPJKM7CzoS/j0NC+W6QkQANgCXichDVfWuyXq6IvI04OWquqWjrPbnjSKykZgwcAo93ctF5BzgfFU9LfrSZkwYKCKLJoliKmPuy9gfs+2nU9z6ZOCyafppSmw5xyT3wjSIyNo4ckD2JPz7paoeoaobVXUjFlXwwGk6MlhmiL1e2URkXcfuASUMBNoJAxMJIPXTpUgyutOzDviBiFwJ/BzzlU37ay8i60Xkgo79Xmzm+LxJl75LRH4Vy30s8JcHUKe/x0YtiUSb1E+XGMm9MA3TJfybdM3Gjs930hFdoKpjwOop7kkJAxNzRuqnS4/UMIlEItFFktFNJBKJLpKMbiKRSHSRZHQTiUSiiySjO8fIHImYx+vOEZFrRGRURG4UkUfOX80TiUQ3SNELc887gBepajMuaPihiPxSVS+F2YuYi8gTMYGbPwYuwUKDEonEEieNdOeYORQxfxvwt6p6saoGVb1DVe+Y+xonEolukka684CIfAh4IdAD/JIpRMzj8szp7s+ABwPfEJEbgAbwNeC1qjo+6doXAy+OuztOPvnkmw+x+muAbamMRVdG4xDvTywSRHU2qm+JAyUazocBj8HcBHXMAD8pCoNsAv5UVf93invXA3cAlwJPAUrg68APVfWNM3znL1T1wYdY71TGYVpGYnGQ3AvzhKp6Vb0QWw75Usxd8FlVnc1ItD2a/VdV3ayq24D3kPR0E4klTzK680+O+XQfD/yFWDr1u9hXxHw3qjqIiZSk15BE4jAjGd05ZBoR82dherhtEfMz4nYn8BJMoHwqPgWcG8tcCbyK/cvrffSQHyKVcTiXkVgEJJ/uHCIiazEt0dOxH7RbgPer6semuHYTHT7dKGL+SFV9ctwvsAwSz8aE0L8MvE5VJ7rwKIlEYr5Q1bQdRhvwSSxty1Udx/4QuBrLbfXgjuMF8BngV8A1wOs7zj0Gy0D8rrj/NOBrHedfD9zQsf8U4BtdqvdzsFRH7S0AZ8xQ9mswV82ajmP/FJ/v0XH/v4Cnd5y/DnhTx/5Xgd9f6P+/aVv6W3IvHH58GvidSceuAn6ffdX8/xCoq+r9sGywL4nZLMAm/x4JZHGRx0VYNEabhwFDInJE3H84lk5o3uutqp/XPVlsnwdsUtXLpyp0qoy402S0vSg+AyKyGksmOvl5LzrAZ0ok9iEZ3cMMVf0xsGPSsWtU9bqpLgf6RCTHYopbwFA85+L5gLmh7gZ2iciJ8fzR2Ojv4XH/4RyCUTrAeneyT9aDSUyVEXeqjLY/Ye9nOR9YK8bxwLhOn3khkZg9Cz3UTtvMG/AK7DW4CXx60rka5kPehBmQx8TjG7FRYh34PyzOVzHf8JM77u/FUrv7eH4TZvjuxkaXVwFXYIa4wgzVz4CTgS9hk4Pvxoxvu4xbgTHgB8BxwP9go8ZmxzWD0zxHO2LjFmAYS/vydCwp59hUZWAJPk+boi2+FuvdjM8yNkU9KuBxsZw6sDM+y854/jvYCsL25wN5ls3x/FBsz//ElnIfSBntZ9mrTeM1D4zPNQJsAV45Qx/6I8x9tLtNF7pf35O3NNJd/NwJvB3zeU7FhViq9qlGYa/EQtP+FPgY9g/4tR3nH4L9g30cZmj6MV/uccANmM+3Xc6zMJ/rqdhI8KeYJsRDMePwm1j+J4BV2A/Ff6jqk1W1P97/B/F5Nk3zHH8bP38QWBbr+gXgWswVsFcZInImMKaqV3WU8Vzsx2VFrNNRWKqZrdiPUWcZv1LV7wOoLd3+NWbgbortcjtm9MeAXx/gs7w7tuk5sT2HgU8dYBkvn6pNRWQN8C3gI1jWhxOxH4V9EJGjgc8B/6+zTTvcQolus9BWP22z2zDD++kZzt/OviPdf2PPRNjbgUHM/9m+54PA8zruvwD4o7j/QMzd8K6O6y/ARmcfAx4aj/0M+4f8HWz01q5DH7bI4z4d92/ERstfmeYZnk7HiD0euxu4jDiR1lkG5jp4w6QyTo7nX4oZ2k1xqzCDetR09cAm116HGbfbsQnG9nP92aRrZ3yWKf6fPBAYPpAysOXdU7XpR7CFNrPpN2cCWycduxt42EL36Xvqlka6hzefAM6Ky4pzbP1+p9/1VuBxYkIQgo1ir43nHoXN4J8Vkxn2AsdiRvGR2JJmsOiBP8cMdNkuWFVHsVf/+3Z83/Mxd8bINPW9Mpb/8Bjn/HTM6HRe31nGH2Jujk7OxH4YHof5pYeBV7N3Rtzp6vETLHb6CszHvQtbUZiz72h0f88ymUdhkRgHUsZ9mbpNHw7sEJGLRGSriPy3iBw7TRm/AK4RkadOatMrZ1nvxByTjO5hhoh8EXv1PxnTa8gxHYe/wozuE0Tk2/HyD2IuhauAI4D/UdUrReT+wN9goVa3xvuHMOGWLcA2VW0bg58CJ8TzYVJ1dgEDHfvPxwzadPX+CWb8/x4z4v8R63wm8M1Y73YZRwG3q+pNk4o6Nz7zKLAe84l/hj2ukpnqcVF8lsuwV/EPYKPlFvuuDpz2WaZ4tnZ7vnbSqf2V0c/UbboWeAHm9jkWuJlpJhNV1QP/jrlpmvHvS6IBTywECz3UvidvwA/ZI/84ebtw0rX7uBf2dz/weSz+dBWm83sZ8LMDKGNrx/31eO3dB/gc18RrHoGN6N45+TlmUcYv56CMm+agjAu7VUa85n3xmsd0HPsVNur+VMex1bHs5VOU8QRsIvLB2CDrIdgk3xkL3f/vqVsa6S4gqvoYVZVptkccyP3YaPSxk+4/HfsHvQN75b4WeGiciNmrjHj/OVjkwEvjsa3t+9Ummn4OrOm8v10G9lreatcBG6WNA8+Il70AOI+O1+UpyngtNgn22I7n+jr24zHbMtrZOJ7QUcb5mAGbbRlXYCPMQWxk3/7837MtI37vXfH+l07x/3XGMiJX0zFCF5E+TMfjCvYeebc/T6UXegbwY1X9hZou888xP/yMWUsS80cyuoscEclFpIHFlmYi0ohxte3z9XgeoBbPt//x/Rx4QZypLoBTsNnynVPc77CJqY/ErX3/80VklYj0Y6OkcWBkch2wSbYC8xEvx16nr1TVa0WkB/O/fn6m58CMSQ24Tzz/AMx/fOUBlPEz7AfmeSLSJyJnYREZ3z6AMs4mGkvMtbFbJ2O2ZYjICZgb4L+BT3f8P+EA6jFlm2I/IM8QkTPicvE3YyPo3f9fO/g58EgROSOWu7tNp7g20Q0Weqidtpk34K3s+4r61o7zm6Y4vzGeW82eyal97sdGUndPcX4Ue60dwQzD6FRlxPufM00dLu6ox7PYE/t6IGXcCrz6EMv4DfCMQyzjdmzkfChljAIjh1hGZ5u+FHs7aY/Aj+noE1cDz+nYfwUWAjiMhcO9eqH79T15S4I3iUQi0UWSeyGRSCS6SDK6iUQi0UWS0U0kEokukoxuIpFIdJFkdBOJRKKLJKObSCQSXSQZ3UQikegiyegmEolEF0lGN5FIJLpIMrqJRCLRRZLRTSQSiS6SjG4ikUh0kWR0E4lEoosko5tIJBJdJJ/ppMiJaprVbf1l6dgm708+xizumWp/ttdMVeFZFru/au7vntmUOZf1OJQyd//V+NmkPCXui0zab5/v+Nv+vOe47i568rnp9ufi77712P/3zkc9Dub596mHtv9acTI5D4R2fJ7u72zOde7PdOxAvm+2ZU5VzwO5fpZlKqBqG+3P+3nU/e1PV/XZ7m+Gb6vq7zAFMxpdGANezp6MIXn83L6tmObYdNfMVRmTEEx/v/OWvONYe799rLPI9vXZfo4daJnTfc90ZeRT7M9UxmzLzGM3yD3kHpdVtlt4stw2gCyvyHNP5uI+Pm4VOfses6/xu4/tOT/19fNTRnu/2l3OzGXM9L3zVcYU13tPVsVjPpBVEP+3IB5LFh/3qbA8GFPtz+aa6fYP5Z6Zypjumpn2D/Ye7HhZQVVBGa+pqngsXlLGrXO/Yk+epGqG87O5Zqrzb7UkrlOS3AuJRCLRRZLRTSQSiS6SjG4ikUh0kWR0E4lEoosko5tIJBJdJBndRCKR6CLJ6CYSiUQXSUY3kUgkukgyuolEItFFktFNJBKJLpKMbiKRSHSRZHQTiUSiiySjm0gkEl0kGd1EIpHoIsnoJhKJRBcRVZ3+pMi3mEEXsousAbYtdCUWCakt9pDaYg+pLfawGNpi23Qi5jMa3cWCiPxCVR+80PVYDKS22ENqiz2kttjDYm+L5F5IJBKJLpKMbiKRSHSRpWJ0P7rQFVhEpLbYQ2qLPaS22MOibosl4dNNJBKJw4WlMtJNJBKJw4JkdBOJRKKLLEqjKyKrROS7IvKb+HflDNdmIvJLETm/m3XsFrNpCxE5RkR+ICLXiMjVIvLKhajrfCAivyMi14nIDSLy11OcFxF5fzx/pYg8cCHq2Q1m0RbPiW1wpYhcJCKnL0Q9u8H+2qLjuoeIiBeRZ3azfjOxKI0u8NfA91T1JOB7cX86Xglc05VaLQyzaYsKeLWqngL8FvByETm1i3WcF0QkAz4IPBk4FXjWFM/1ZOCkuL0Y+LeuVrJLzLItbgYerar3B/6ORT6hdLDMsi3a170T+HZ3azgzi9XoPg34TPz8GeDpU10kIhuA3wU+3p1qLQj7bQtV3ayql8XPw9iP0NHdquA88lDgBlW9SVVbwJew9ujkacC/q3ExsEJE1nW7ol1gv22hqhep6mDcvRjY0OU6dovZ9AuAc4GvAlu7Wbn9sViN7pGquhnMoABHTHPdvwCvA0KX6rUQzLYtABCRjcADgJ/Nf9XmnaOB2zr2b2ffH5PZXHM4cKDP+SLgf+a1RgvHfttCRI4GngF8uIv1mhX5Qn2xiPwvcNQUp944y/t/D9iqqpeKyGPmsGpd51DboqOcfuyX/VWqOjQXdVtgZIpjk2McZ3PN4cCsn1NEHosZ3UfMa40Wjtm0xb8Af6WqXmSqyxeOBTO6qvqE6c6JyBYRWaeqm+Or4lSvB2cBTxWRs4EGsExEPqeqz52nKs8bc9AWiEiBGdzPq+p581TVbnM7cEzH/gbgzoO45nBgVs8pIvfH3G1PVtXtXapbt5lNWzwY+FI0uGuAs0WkUtWvdaWGM7BY3QvfAF4QP78A+PrkC1T19aq6QVU3AucA31+KBncW7LctxHrWJ4BrVPU9XazbfPNz4CQROV5Eatj/529MuuYbwPNjFMNvAbva7pjDjP22hYgcC5wHPE9Vr1+AOnaL/baFqh6vqhujffgK8LLFYHBh8RrdfwSeKCK/AZ4Y9xGR9SJywYLWrPvMpi3OAp4HPE5ELo/b2QtT3blDVSvgFdjs8zXAl1X1ahH5cxH583jZBcBNwA3Ax4CXLUhl55lZtsXfAKuBD8U+8IsFqu68Msu2WLSkZcCJRCLRRRbrSDeRSCQOS5LRTSQSiS6SjG4ikUh0kWR0E4lEoosko5tIJBJdJBndRCKR6CLJ6CYSiUQX+f/FQG9MANZFEQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check data read in if flag set\n",
    "if input_check:\n",
    "    sar.plot(faults=None,figure=127,norm=[-0.5,0.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eligible-necessity",
   "metadata": {},
   "source": [
    "Now we can estimate the covariance of the dataset using the semivariogram method. First we create a image covariance object from the InSAR object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "excellent-birth",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------\n",
      "---------------------------------\n",
      "Initialize InSAR covariance tools Covariance estimator\n"
     ]
    }
   ],
   "source": [
    "# Covariance first estimate on original data\n",
    "covar = imcov('Covariance estimator',sar, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "computational-hamburg",
   "metadata": {},
   "source": [
    "We will mask out the area close to the main rupture to look at the rest of the scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "whole-sitting",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Masking data set Covariance estimator\n",
      "     Mask: -118.0 <= Lon <= -117.25 || 35.55 <= Lat <= 35.95\n"
     ]
    }
   ],
   "source": [
    "# mask out high deformation above earthquake rupture\n",
    "covar.maskOut([[-118.0, -117.25, 35.55, 35.95]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "strategic-sensitivity",
   "metadata": {},
   "source": [
    "We use the `computeCovariance` method to sample the dataset with a random set on 0.002 of the pixels, estimate and remove a ramp, calculate the semivariogram at distances of every 2 km out to 100 km, convert the semivariogram to covariance, and estimate a fit of an exponential function vs. distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "certified-quantum",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing covariograms\n",
      "Computing 1-D empirical semivariogram function for data set Covariance estimator\n",
      "Selecting 16028 random samples to estimate the covariance function\n",
      "Estimated Orbital Plane: -3.4665518852367864e-05xy + 0.13375223760488297x + 0.01550056454564969y + -59.674849380877234\n",
      "Build the permutations\n",
      "Digitize the histogram\n",
      "Fitting Covariance functions\n",
      "Dataset Covariance estimator:\n",
      "A prior values: Sill | Sigma | Lambda\n",
      "                 0.000850 | 0.004217 | 39.788291\n",
      "Dataset Covariance estimator:\n",
      "     Sill   :  0.00039557243022636426\n",
      "     Sigma  :  0.07061470016722794\n",
      "     Lambda :  39.78829171795961\n"
     ]
    }
   ],
   "source": [
    "covar.computeCovariance(function='exp', frac=0.002, every=2.0, distmax=100., rampEst=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bridal-globe",
   "metadata": {},
   "source": [
    "Now we can plot the results with the `plot` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "stylish-circulation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEICAYAAACj2qi6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3yElEQVR4nO3dd3hUVfrA8e+bSQEySu9IU7CAkgASxRWQIiQGUVdWQda6C66Crq76gwUbsLZ1BVexYFcQVhEpkVAUASmBBBKKKL33DpME0s7vjzuBIcxMJoVMMvN+nuc+mXvn3jvn3sC8Oec951wxxqCUUkr5IsTfBVBKKVVxaNBQSinlMw0aSimlfKZBQymllM80aCillPKZBg2llFI+06ChgpKIfCAiz/u7HBeDiPwqIl38XQ4VmETHaaiyJCL9gaeBq4BTQBrwL2PMYn+Wq6ISkc+B3caYEWXwWV2ACcaYRhf7s1T5pTUNVWZE5GlgLPAKUBdoDLwH9CnjctjK8vOURURC/V0GVQqMMbroctEXoCrgAPp62ScCK6jsdS5jgQjne78B8S77hgKHgbbO9W+B/cAJYBHQymXfz4H3gVlAOtDduW208/3qQAJwCDjmfN3I5fgFwChgCVbtaC5Qy+X9PwBLgePALuBBl+t5E9gJHAA+ACp7uf6Hndd5DJgDNHFuF2AMcNB5fWuA1sBAIBvIct7bmc79twPdna9fct6bCc6yrwVaAsOc59sF3OpShoecZTgFbAUGObdHAplAnvOzHECDQn5nXYDdwP85fzdf+fvfoS4lX7SmocrKjUAl4Hsv+wwHbgCigDZAByC/2WUS0M9l357AYWPMKud6ItACqAOsAiYWOHd/4F/AJUDBprAQ4DOgCVbtJxN4183xDznPHw48AyAijZ2f/Q5Q21n2NOcxr2N9QUcBVwANgRfcXbiI3AH8E7jLeZ5fnNcMcCvQyXmuasA9wBFjzHjndb5hjLEbY3q7OzfQG/gKKzimYgWkEGd5RgIfuux7EIgHLnVe7xgRaWuMSQdigb3Oz7IbY/bi/XcGUA+ogXVvB3oon6pI/B21dAmOBbgP2F/IPluAOJf1nsB25+srsP76reJcnwi84OE81QADVHWufw58WWCfz3HWNNwcHwUcc1lfAIxwWX8MmO18PQz43s05BKtWc7nLthuBbR4+MxF4xGU9BMjA+rLtCmzE+nIOKew6uLCmMc/lvd5YtQSbc/0S572q5qFc04Anna+7YOVPfP2ddcGqBVXy978/XUpv0ZqGKitHgFqFtGs3AHa4rO9wbsMYsxmr2aS3iFQBbge+BitHISKvicgWETmJ9aUJUMvlXLs8faiIVBGRD0Vkh/P4RUC1ArmP/S6vMwC78/VlWF+cBdUGqgArReS4iBwHZju3u9MEeNtl36NYgaehMWY+Vs1nHHBARMaLyKWerseNAy6vM7FqaLku6+Rfj4jEikiSiBx1liOO8+9jQR5/Z06HjDGni1BWVc5p0FBlZRlwGrjDyz57sb488zV2bsuX30TVB1jvDCRgNR31wcpVVAWaOreLy7Heugn+A7gSiDHGXIrVFFTweE92AZe72X4Y6wu5lTGmmnOpaoyxu9k3/zyDXPatZoypbIxZCmCM+a8xph3QCquZ6lkfrqtIRCQC+A4rD1PXGFMNKw+Ufx/cfVZhvzPtnhlgNGioMmGMOYHVnj9ORO5w/nUf5vzL9g3nbpOAESJSW0RqOfef4HKayVjt+3/DWctwugQ4g1WbqYLVO6soLsH6gj8uIjWAF4tw7ESgu4j8SURCRaSmiEQZY/KAj7ByAnUARKShiPT0cJ4PgGEi0sq5b1UR6et8fb2IxIhIGFaT12kgv6ZwAGhetMv1KBwrsX0IyBGRWKz7ne8AUFNEqrpsK+x3pgKMBg1VZowxb2GN0RiB9cW0CxiM1W4OMBpIweodtBYroT3a5fh9WDWWjsD/XE79JVazyB5gPZBUxKKNBSpj1Q6SsJqRfL2mnVhNOP/AalJKw0oIg9VraDOQ5Gz2+hGrRuPuPN9jJc4nO/ddh5V4Bisp/RFWr6odWMHxTed7nwDXOJu1pvlabg9lOAU8AXzj/Kz+wAyX93/HChJbnZ/XgEJ+Zyrw6OA+pZRSPtOahlJKKZ9p0FBKKeUzDRpKKaV8pkFDKaWUzwJqArFatWqZpk2b+rsYSilVoaxcufKwMcbTwNPzBFTQaNq0KSkpKf4uhlJKVSgisqPwvSzaPKWUUspnGjSUUkr5TIOGUkopn2nQUEop5TMNGkoppXwWUL2niis3N5fExERSU1OJjo4mNjYWm00fI62UUgUFfdDIzc2lZ8+eLF++nPT0dCIjI4mJiWHOnDkaOJRSqoCgb55KTExk+fLlOBwOjDE4HA6WL19OYmKiv4umlFLlTtAHjdTUVHIcDsZy7pmW6enppKWl+a9QSilVTgV90IiOjqZT5coMwnp6TAwQGRlJVFSUfwumlFLlUNAHjdjYWHI7dqR75cpkA4uAfzVoQGyvXv4umlJKlTtBnwi32WzMmTOHxMREpi5bxp/nzuWJlBR44AFy33+fxEWLtFeVUko5BdTjXtu3b29KPGFhXh68+irm+efZXqUKfzSGtMxM7VWllApYIrLSGNPel32DvnnqAiEhMHw4y19+mUvS01mQkcEftVeVUkoBGjQ8mge0BdYD3wLvAbkOh/aqUkoFtaDPaXgSHR3NMbudTg4H/wKeBW4OCeFwzZo6glwpFbQ0p+FBwZHid0RE8GluLlXDwvhPkya8vHMn6RkZmutQSlV4RclpaE3DA9deVWlpaURFRXFJmzYciY/nmTVrqA/8DTjlkuuIj4/3d7GVUuqi0qDhhc1mIz4+/rxg8MFdd3FmzRpewhoIeB+Q7BxBrkFDKRXoNBFeRFHt2jHWbqczEA4sAUaFhRF97bV+LplSSl18pRI0RKSXiGwQkc0iMtTN+yIi/3W+v0ZE2hbh2GdExIhIrYLv+UNsbCwxMTGstttpA0wJDWV4VhZxr7wCmzf7u3hKKXVRlbh5SkRswDigB7AbSBaRGcaY9S67xQItnEsM8D4QU9ixInKZ872dJS1naSmY67BHRZF36hQhjz0GUVHk/ec/zGrQgNS0NO1ZpZQKOKWR0+gAbDbGbAUQkclAH6whDvn6AF8aq6tWkohUE5H6QNNCjh0DPAdML4Vylhp3uQ5uvhlz//2EPPooYrPxbm4uGXa79qxSSgWU0mieagjsclnf7dzmyz4ejxWR24E9xpjV3j5cRAaKSIqIpBw6dKh4V1AaGjXih7//naHh4XTLzWUd0EtHkSulAkxpBA1xs63g4A9P+7jdLiJVgOHAC4V9uDFmvDGmvTGmfe3atQst7MWUuno1b2Rn0w7YjjWS/FOHg42LF/u1XEopVVpKI2jsBi5zWW8E7PVxH0/bLweaAatFZLtz+yoRqVcK5b1ooqOjiYyMZD1wIzAMuB14/P334ZtvIIAGUiqlglNpBI1koIWINBORcOBeYEaBfWYA9zt7Ud0AnDDG7PN0rDFmrTGmjjGmqTGmKVZwaWuM2V8K5b1o8ntW2e128kR4127nsZgYwlu2hHvugb594cABfxdTKaWKrcSJcGNMjogMBuYANuBTY8yvIvKo8/0PgFlAHLAZyAAe8nZsScvkL+5GkcfGxiLGkPfvf2NeeIEziYlsGDSI6954A1uojq1USlUsOvdUGcifx+r4smW8k5HBjcDKGjWIWrYMW8uW/i6eUirI6fM0ypnExESWL1/OyowMbgIeA1ocPQrXXguvvw7Z2f4uolJK+USDRhlITU0lPT0dsLqMvQ9cA2y6/HIYOhTat4cVK/xZRKWU8okGjTKQ36vK1Qm7nc1vvAFTp8Lhw3DDDTBkCBw/7p9CKqWUDzRolAHXXlUigt05Ujw2NhbuvJPctWvZHhdH3rhxnGnWjLzPP9fuuUqpckkT4WUk/2l/rr2qbDbbeQ97aulw8EFICNfn5WFuugl57z247jp/F10pFeCKkgjXoOFnCQkJ9OvXD4fDAVhD5P8WEcGY8HDCMzJg8GB4+WWoWtW/BVVKBSztPVWBuCbJwZkoz8riv4MHw1//Cv/9L1x5JXz6KeTl+a+gSimFBg2/c5ckj4yM5KqOHeH9961eVc2awSOPwPXXwy+/+KmkSimlQcPvvCbJweqOu3QpTJwIBw9Cp07wpz/B9u1+LbdSKjjpPBZ+5mnqkfwkeWJiIqmpqdYDndavx/bWW9aAwBkz4B//gGHDwG7392UopYKEJsLLKddeVenp6URGRp57oNO+fdagwIkToW5dePFF+MtfICzM38VWSlVAmggPAPlTjzgcDowxOFwf6NSoEUyYAElJVpL8scegdWtroGAA/RGglCp/NGiUUwV7VQGkp6eTlpZ2bkNMDCxYYDVVhYbCH/8IN90E+tAnpdRFokGjnPLUqyoqKgqwmq8SEhIYNXo0CSLkrloFH38MO3bAzTdDnz6wdq0fSq6UCmSaCC+n8ntVFcxpxMbGes939OsHY8dayfI2beDee63BgS1a+PuSlFIBQBPh5ZinqUcKjiIHsNvtTJo0ifj4eGvD0aPw739bgwPPnIEHHoAXXoAmTfx0NUqp8koT4QHCZrMRHx/PiBEjiI+Px2azAT7mO2rUgFdfha1brdlzJ060ahuDB8Pego9wV0op32jQqIC85TvO5jpGjSIhIYHcWrVgzBjYtAkefhg+/BCaN7cCye7dfroCpVRFpc1TFZCnnMasWbOIi4tzn+tw1lLYts2qgXz2GYSEWIFk6FBttlIqiOkst0HAXb4jMTGx8FxHvh074LXX4JNPrLEdDz5ojS5v3rxsL0Qp5Xea0wgC7vIdheU6zmu6WruW3HffhS1b4NFH4auvoGVLuO8+WLPGD1eklKoINGgEkMJyHT179qRfv368+OKL9OvXj549e5LboAG8846VMH/qKWugYJs2cNtt1oy6AVQTVUqVnAaNAOJtxlyv05IANGhgddHduRNGj4bkZGtG3ZtusgKJPstDKYUGjYCSP2PupEmTGDlyJJMmTTqbBPfWdHVes9WSJeQOHWpNvf7uu7BvnzW6vFUrGD8eMjP9c3FKqXJBE+FBwtOAwAkTJvDOO+947HGVe+YMq4cPp/7XX1N/3z5MzZrIY49ZkyTWq+fHK1JKlRZNhKsLeGq6Ajw2W+Xm5tLzttvo/OGHNNy3j16VK7NUBDN6NKZJE3b26MGHgwdb40Fyc/18hUqpsqBzTwUJTw97euWVV7z2uMoPKABzMjNZYrPx/bhxOF55hR4//sigH3/kl/ffZ3SrVoxITsYWEXHhw6M8PVTKuV0pVXFo0Agi+d10Xcds5Pe4cm22yu9x5SkP8sXSpUw7fpww4BHgsbw8Xly7lsxGjYh44gnumTePOc5jizTwUClV7mnzVJDz1uPKUxdeYwzp6ekcA94ErgD6APuqViXkhReY+MsvjHM4iHFp7ho9erTX3lsXTH+izV1KlUta0why3p5R7ml69r59+zJ9+vSztZM8YL7dzvqxY5k9dy5577zDA8D9wBpgvMNB2sKFHpvBYmNjPU71DmiTllLliPaeUl65m64E8Pglnz+VCQ4H/YBBQDsgKyyMycbwbk4Oyc5z509xAhSvZ5fmSJQqFUXpPYUxJmCWdu3aGVU2cnJyzMyZM82oUaPMzJkzTU5Oztnt3bp1M3a73YiIsdvt5m/XX29yHn7YZNhsxoBZBeYf4eHmjk6dTE5Ojhk5cqQREQOcXUTEDBgwwNjt9vO22+32s59X8HO6det2thxKKd8BKcbH71ltnlLF4i6pnr/dU3NX+Jtvsnb4cOpNn86be/dikpKQ/v25tVUr3qxShZMuzVeuuRNXnnp2ueZI8kfAaw1EqdKnzVPKP9LSrOnZJ0yAo0c5GBHBZ3l5fJydzX5nMn7IkCEMGDDA7ay9qampvPjii7j++xURXnrpJRYtWqRNWkoVQVGap7SmofwjKgrefhveeANmzKD2J5/w3Ny5/B9wtFEjqt15J6ZjR4/PSQfcdhXOycnxWgPRhLtSJVMqQUNEegFvAzbgY2PMawXeF+f7cUAG8KAxZpW3Y0Xk30BvIAvYAjxkjDleGuVV5UhEBPTti/Ttaz1J8OuvqfHVV9ZjaZ96inlxcax84gnmhYVxbfv2hfbsstlsRW7SSkhI8JpwV0q58DX54WnB+rLfAjQHwoHVwDUF9okDEgEBbgCWF3YscCsQ6nz9OvB6YWXRRHiAyMszJjXVmKeeMqZuXWPAmGrVjHnkEWN+/NEYl6R7wWT8zJkzPSbPi5NwVyoYUIREeGkM7usAbDbGbDXGZAGTscZ6ueoDfOksXxJQTUTqezvWGDPXGJPjPD4JaFQKZVUVgYjVfPXWW1btIzER4uPhf/+D7t2hYUN44glsK1YQf9tt5z2IqiSDFV15fHiVDjxUQa40mqcaArtc1ncDMT7s09DHYwEeBv7n7sNFZCAwEKBx48ZFKbeqCEJDoVcva8nMhB9+gEmTrGna33nHerZ5377Wcv31pTJYES58eJU2XSllKY2gIW62FeyS5WmfQo8VkeFADjDR3YcbY8YD48HqPVVYYVUFVrky3H23tZw8CdOmweTJVkL9zTehcWP44x+x9e1LfFycz92BAY8Jd9eHV4F27VWqNILGbuAyl/VGwF4f9wn3dqyIPADEA92c7W5KWS69FO6/31qOHYOZM2HKFBg3DsaMsZqw7roL7rwTbr7ZqrHgeXyJp9qJp0kbV61axdixY7UGooKPr8kPTwtW4NkKNONcMrtVgX1u4/xE+IrCjgV6AeuB2r6WRRPhypw4YcyECcb06WNMpUpWEr1GDWPuv9+YqVONcTiKdDpPifXnn39ek+cqYFCWiXBjJasHA3OA34BvjDG/isijIvKoc7dZzuCwGfgIeMzbsc5j3gUuAeaJSJqIfFDSsqogcOmlcN99VtPV4cMwdaqVRJ8506p51KplPb7244+tR9kWwlNi3VvXXk2cq0CmI8JVcMjOhsWLrWAybRrs3Gltb9/eCiq9e0N0tNVzqwB3kzbmT8xY1EkWlSqPijIiXIOGCj7GwLp1kJBg1UCSkqxtDRpAXJzVU6t7d6ha1eMpPPWq8jb1ScE8ilLlhU4jopQ3InDttdYybBgcOmSNBZk5E7791mq6stmgY0eIjbWCSJs2EHKuNbc4j8+Nj4/Xua9Uhac1DaVcZWdbNY/Zs61Akppqba9bF3r0sJbu3a1aiRsJCQlum60mTZrkde4rDRzKn7R5SqnSsn8/zJljBZEff7SS6wCtWp0LIp06gd0OeG62cn1AlaeAojUQ5S/aPKVUaalXDx54wFry8mD1apg3z1refx/GjrXGgMTEwC23YOvalTnTp5P488865kMFpNKYe0qp4BASYvWweu45K2gcOwZz58Izz1jNWq+8Al27YqtZk/i33mKEMcRfeim27GwAj3NfuU7nbow5b9S5UuWN1jSUKq7Klc81UQGcOAG//ALz58PPP8MLL1jbIyIgJoa4P/yBv11xBV9s2sShjAyfpnPX5LkqbzRoKFVaqla1xnzkd609etQaG7JwISxaRMhrr/FGXh6v2Wzsa9iQM+3a0WTAAH7KzHT7QCmdMFGVR5oIV6qsnDwJy5ZZQWTxYlixAs6cAWBP5cosyM5mYU4OaZUrU+3GG0mcO1eT56pMaCJcqfLo0kuhZ09rAStgrFoFS5bQYPFi7l6wgPtOnIDMTExyMtKzJ1VzcujscJAEHHGeRpPnyp+0pqFUeWEMbN5sjRNJSoJly8hbvZqQvDzAmrhtBbA6PJwW/fvzz2+/5ZBLLsR15LnmQVRRaE1DqYpIBFq0sJY//xkAc/IkT99yC5esW8d1WVl0EqF/VhZ8/jkPAmuBZGAlsMrhYG1Kig4iVBeV1jSUKucumDAxKopVH37Iz6+/TpvsbK4Hajj3zQsN5VTjxny/cydJOTmswgosoZoHUV7oiHClAtx5vaocDq6pUoW+zZvzfFwc26dModrWrWcDSS7wO5B37bUszcwkYfdukk6f5rRzmvf8Gog2aQUvDRpKBQF3U7bbbDZr/qt776VWejpRQDTQ3mbjpshIqp48efb4PcCvNhvN77iD5n36MHDcOL5bt44TLmNI5syZA+AxmGigCQwaNJQKYp7Gdtx88828+9JLXIcVSK4FrgOutdkIdT4oKgurVrIO2BQeTq9nnuHdBQuYvno1DjfBxFPuBDwHGlX+aNBQKsgV5cFRk7/6il0//sgv48ZxHdAKaA00dTlfOtazl9cDW8PDiX/uOU5ddhl3Pv00Jwv04CrsQVRaOyl/NGgopS5Q1Bl460dG8tc//IE9c+acDSRXA41cznka2ID1rOYNzqV5XBwfLVzIQTfdgbVnV/mkQUMp5ZanPEhRnkTYIDKS/738MvZdu1j43ntckZ3NlUAzwPVrfydWU9cm5xLz5z9T9+abueuppzjhIaBoDcQ/NGgopYrMXUAB73kL1/dqVKnCHa1bM7BzZxLHjqV5VhZXAi2Bai6fkwVsxQokW5zLDQMGMGfLFmauWXNBIl6btC4+DRpKqVLjqXbi6T0oEGiqVOHWtm35ZvRo1k6dyk/vvUfT7GwuB64A7C6flQPswAoqu0JDufG++2gZG8vg//yHhPXr2VOEnl3Kdxo0lFJ+5XMzWJUq9IqOJrZFC3757DOaYwWS5s6ldoHzHge2A7tsNq6JiyNxwwZ+2bmT30+f5nBkJFffcIPWTopBg4ZSqtwqSs+ubz/5hJ0LFzL7vfdojtWjqylW/uQKm40IZ1fhfMeBkKZNuaR1a6anpZF88CBbsrI4VLkytdu2ZeL8+WCz6biTAjRoKKUqlOI8W73P7bcz9+uvzwaSJs7l1iuvpN6ZM+Rt335eLgUgLySEQ+HhbM3KYnteHgfDwght2pRH//UvaNiQ+557jtlpaZwMsgGOGjSUUhVOafTsyu+JlZqayosvvsglxnAZcBlWQIlt3Zr033+nfk4OlwENgcoFywHsB/YCB2w2WvXowc8bN7Jizx62nznDscqVadi+Pd/Onw8iRR7gWB6DjAYNpVRAKWrPLo+1kz59+Prrr3H93qsJvPLYY1xy4gQLJ06kIdAAK6A0BJqGh3NJVtYFZcoLDeVMtWqsO3qU3Xl57McKNsciIrjv6af5at485q9fz7aMDEKd83zNmjWLuLi4cjfwUYOGUioolGbtBPAYaKZOnEh9OLs0AP50883YT53iQFoa9YF6XJi4z3cCOCRCWKNGpO3dy57cXA4AB4GTlSrx+IsvEtO7N3c/9hg/rVzpdsqWixlMNGgopYJeScedFCfQhALNqlThT506kTZ7NvWAukAd58+rqlUj7Phx6gK1PJQ7BzgMHAKO2mxcERNDyo4d/HrwIPuyszkVEUGda67h1fHjsdWtS261aiQuXFiigKJBQymlPCjxuJMSBJqnnnqKMWPG4HA4CMVqGmtauTJj//lP9q1ezcIpU6iDFVBqO5fml15KpZMnz051704GMAn4e4Hp7n2lQUMppUpRaQUabzmNwvIwIcZQEyvQ1AKG9OtHy5o1mTp+PPasLNYDn3P+Y399pUFDKaX8zFu+pbR7ibl+j4sII0eOZMSIET6XVZ8RrpRSfmaz2YiPj7/gL35v2/NrHAVrLTExMRcEk/z3IiMjzwsokZGRREVFXbTr0pqGUkqVc0WtnWhOw0caNJRSwcZbvsVXGjSUUkr5rChBI6SUPrCXiGwQkc0iMtTN+yIi/3W+v0ZE2hZ2rIjUEJF5IrLJ+bN6aZRVKaUCSW4uJCTAqFHWzwJzOJa6EifCRcQGjAN6ALuBZBGZYYxZ77JbLNDCucQA7wMxhRw7FPjJGPOaM5gMBf6vpOVVSqmKJjcXEhMhNRWioyE2Fmw2a3vPnrB8OaSnQ2QkxMTAnDnW+xdDafSe6gBsNsZsBRCRyUAfrGfQ5+sDfGmstrAkEakmIvWxJqf0dGwfoIvz+C+ABWjQUEoFMHfBATwHhsREa3t+5ymHw1pPTIQiDNMoktIIGg2BXS7ru7FqE4Xt09DD9vxj6xpj9gEYY/aJSB13Hy4iA4GBAI0bNy7mJSilVOnyVjsoSq1hyJCCgSGHpKRjjB9/mJUrj+BwHAaOAI2AnqSnQ1pa+Q4a4mZbwey6p318OdYrY8x4YDxYifCiHKuUUr7w9EXv6T1wHwBmzYK4uPO3X3/9GT777BDff3+QJUsOcfr0QeAQDschFi48xIYNh3E4DmHNRnUYOEZ6Ojz2WMFS/hHoSWQkXMRhGqUSNHZjTVefrxHWVPS+7BPu5dgDIlLfWcuojzUhpFJKlUhpBQDn1FNu3/vrX9NZtmw/GRnWhOkOx0EWLTpAx44HSEs7SG6uNcetw3GQn38+QdOm7koaSk5ObRyOWthstcnNbYs1gUgtwsNr8eSTtejatSYvvliLdetqkZFRC7vd+vz8cl8MJe5yKyKhwEagG7AHSAb6G2N+ddnnNmAwEIfV/PRfY0wHb8eKyL+BIy6J8BrGmOe8lUW73Cql8pVWABgyBAYMAIfDYDUD7aNSpX0MGrSXY8f2MWnSPrKz98HZJ2rsBxzuikRYWC2ys/PnvM2f/7YOt99em9at6zBmTB0yM/OnKqyK3S5MmADvvOM52Z1/nWlpVg3DNQj6qszHaYhIHDAWsAGfGmP+JSKPAhhjPhARAd4FemFNyPiQMSbF07HO7TWBb4DGwE6grzHmqLdyaNBQKrgUJz9gBYBz57Db4YsvTnPkyB6efHI3mZm7sRo89mCz7aF69b0cPrwH2Adc+DAmuJRzT9qoB9QjOroe69fX48yZ/MnR6xIZWZunnw5lzJgLP3/SJKvs3gJaSQODNzq4TykVUIrTq6hfv/wv5yxgD5Ur7yQ6ehdLl+7C+jt0t8ty2M2n2oGG1K3bgCNHGpKT05D8RzBVrlyft9+uT/Xq9XnooSoXBAFPtQN3OY3SrjUUhwYNpVTAKLzWcBLYDuwgPHwn8fE72LhxB+vW7cAKDvu5sH9NTawU6mWEhjbi3nsbUbNmQz78sBGnTzd0vneJ1wDgrUnLW+3AX4HBGw0aSqkKx1NT03ffZXL//dvJyNgKbAO2Y7Nto2rV7Rw9uh04v9U6NDSCOnUac+BAE3JzG2O1cDemUqXL+M9/LuPbby8jJaVKqQQA13KXpyBQVBo0lFLlVsHg0KuX4ciRI8THb2bt2i2cObOF0NAtREZuJTJyK3v3FuyMWQloSoMGTTl4sBk5OU2xxgk3oUqVJkyaVIfbbgsJ2gBQHBo0lFJ+5a7WcPz4EX77bSODBm1k06ZNZGdvJiTEWnJyTrgcLUBDQkKac8stzalXrzlTpjTnzJlmQHOgbpn1KgoW+hAmpVSpKcq4hq5dz7B582YGDPidDRt+JytrIyEh1pKT49qMZAOakZd3OTbbDbRrdwUrV14BXIFVa6iEMdClCwwbBvv3Xxgc4uOtxVNgsNnO7aNKjwYNpRRQtB5K//vfcX777TcGDlzP5s2/kZ39OyK/Y8w2IM/lrI3Iy7uSkJA/8cgjLTl9uiUTJ7YAmgFhAOTkwNVXw4YN53dFzR/ZbLOd6xHlLjhoYChbGjSUCiJFHdfw0EPHWbr0VzIz1wG/4nCsZ/789dSqtc/lrJWAKzGmHWFh99Ghw1UsWXIV0BKIPPu5TZtaX/jTp18YHO6+G/btu/Dz8wOX1hrKDw0aSgUJb9NoT59+mqVL15OZuQZYh8Oxjvnz1/HTT3tczmAHrsGYW7n11muoXfsaJk68BmiC1dxk1RqaNYPVq93XGmJjrc8salOTKj80aCgVgNzVKBITISnJkJ6+F0jD4VjDwoVraNZsNXv2bCQvL//pPZWAqzGmK9HRrfntt9acPt0Kq+uqYLdbYySg6LUGbWqq+LT3lFIBJr9GkZSUS3r6RiIiUqlbN43Q0FS2bk3j/NHPTbnqquto0+Y6pk+/jtOnr8NKRtsuysA2VT5pl1ulgoBrbeK663Jo2vR3UlNT+O67lcyatZK8vDQg07l3OPXqtebw4WhycqKAKOBa7PaqJZr3SLu1BgYNGkoFsLy8PDZs2Mzdd69g06ZksrOTgTTyA0R4eCRZWdFAO6AtEA1cxcsvh7FokY5rUBfScRpKBZC9ew/w3ntJ/PLLCk6eXMG2bcmcOJE/GK4KVmAYREREO8aMaUeDBi0ZMMB2wSR6bdvC8OE6rkGVjAYNpcqRrKwsVq9eTVJSEsuWLSMpKYlt27Y53w0lJORa6te/l86dr2fGjA7A1eT/N87KgiNHYOBA9z2U8gOEBgZVEho0lPKjo0ePsnTpUpYsWcKSJUtITk7m9OnTADRs2JCmTW9kz57Hycq6AWhLXl5lTpyANm1g/vziDYZTqiQ0aChVhnbs2MHChQtZvHgxS5YsYf369QCEhoZy+eXtiI7+G507d+TRR2+gSZNGjBoFS5eef470dAgN9VybAK1RqItHg4ZSF4kxhk2bNrFw4UIWLVrEokWL2LlzJwBVq1alY8eO3Hfffdxww02MGnU9KSlV2LgR1q6F5GSrthAdbQWEgjWKwvITSl0s2ntKqVJijGHLli3Mnz+fn3/+mQULFrB//34A6tSpQ6dOnfjDHzoREtKZo0db065dyNlBd+eeMmfx5RGgGiBUadHeU0qVkV27djF//vyzgWLXrl0A1K9fn65du9K5c2c6d+5My5YtycsTtwHg5putdVfp6VYNIj5e8xOqfNGgoVQRnDx5kgULFjBv3jzmzZvHhg0bAKhVqxZdunRh2LBh3HLLLVx55ZWIyHnHJiZaASO/RuFwWOsdO7pvgoqKsl5rfkKVJxo0lPIiNzeXlJQUZs+ezbx580hKSiI3N5fKlSvTuXNnBg4cSLdu3bj22msJCQlxHgM//HDhTLKpqe5rFIUltZUqTzRoKFXAgQMHmDNnDrNnz2bu3LkcOXIEEaFdu3Y899xz9OjRg44dOxIREXHBsd5mktWktgoEGjRU0MvLyyM5OZmEhARmzZrFqlWrACt5fdttt9GrVy969OhBrVq1Cj2XpyaoxETP04LroDtVkWjQUEHp1KlTzJs3j4SEBH744QcOHjxISEgIN954I6NHjyY2NpaoqKizTU7uuJt+3FMTlCa1VaDQoKGCxp49e5gxYwbTpk1jwYIFZGVlUa1aNXr16kXv3r3p1asXNWrU8OlcnpqhhgzRpLYKbBo0VED77bffmDZtGtOmTWPFihUAtGjRgiFDhtC7d286duxIWFhYkc/rqRlqyBBNaqvApkFDBRRjDKtWrWLKlCl8//33Z7vEdujQgVdeeYU77riDq6666oLusN4UpRlq7VptglKBTYOGqvCMMSQnJzNlyhSmTJnCtm3bsNlsdO3alSeffJLbb7+dhg0bFuvcxWmG0iYoFcg0aKgKKT9QTJ48me+++46dO3cSFhZG9+7def7557n99tupWbOmz+dzV5uw2bQZSqmCNGioCmXdunVMmjSJyZMns3XrVsLDw7n11lsZNWoUvXv3pnr16l6PdxccwPPYCm2GUup8GjRUubdt2za+/vprJk+ezLp16wgJCaFbt26MGDGCO++8k2rVqp23v6dag7emJk9jKzwNyNNmKBWsNGiocunYsWN8++23fPXVVyxevBiAm266iXfffZe7776bunXruj3O24hsT01N9ep5HlsxbJg2QynlSoOGKjeysrKYPXs2X331FTNmzCArK4urr76aV155hf79+9OkSZPz9ndXo/A2IttTU5OI99qENkMpdY4GDeV3a9eu5dNPP2XChAkcPnyY2rVr8+ijj/LnP/+ZqKh2zJ4tfPmlb01N3qYZ99TUdPfdsG+fPgVPKV9o0FB+cfz4cSZNmsSnn35KSkoKYWFh3H777Tz44IP07NmTsLCwYjU1eZtm3NPcT/kBQWsTShWuREFDRGoA/wOaAtuBPxljjrnZrxfwNmADPjbGvObteBHpAbwGhANZwLPGmPklKavyP2MMCxYs4OOPP2bq1KmcPn2a6667jrfffpt77ulPcnItUlOtfYvb1ORtmvHCmpq0NqFU4Upa0xgK/GSMeU1EhjrX/891BxGxAeOAHsBuIFlEZhhj1ns5/jDQ2xizV0RaA3OA4o3OUn53+PBhvvjiC8aPH8/GjRupVKkaXbs+zEsvPUL79tHFeqJdcacZ16YmpUrIGFPsBdgA1He+rg9scLPPjcAcl/VhwLAiHC/AESCisPK0a9fOqLKRk2PMzJnGjBxp/czJOf+9GTPyzMMPLzSdO/c34eHhBjBVq3Y0ERFfGMgwdrsx3bqdO4/dbgycW+x2Y55/3v32/M/r1s1aFzHnnU8pVTRAivHxe7+kNY26xph9zuCzT0TquNmnIbDLZX03EFOE4/8IpBpjzpSwrKoYijoY7tSpU8TEfMmmTeMw5jegKo0aDeTZZwcxfHhrzjh/ixe7qUkpdXEUGjRE5Eegnpu3hvv4Ge5mhjM+HSjSCngduNXLPgOBgQCNGzf2sUjKF0UZDLds2QbuuGMcP/30OZmZp4D2wKfAPRw/XoUVK7SpSalAUGjQMMZ09/SeiBwQkfrOWkJ94KCb3XYDl7msNwL2Ol97PF5EGgHfA/cbY7Z4Kd94YDxA+/btfQpG6nxFnXfp3GC4XGAW8A4ZGfNITAyjdet7WL16MOcqk4WPhdAn2ilVcZS0eWoG8ABWT6cHgOlu9kkGWohIM2APcC/Q39vxIlIN+AEr97GkhGVUXnjr1uqp2Sgnx0FY2OdkZY0FtgANCQ8fzccf/4Xq1evSr1/RxkJoU5NSFYdYOZBiHixSE/gGaAzsBPoaY46KSAOsrrVxzv3igLFYXW4/Ncb8q5DjR2AlzDe5fNytxhh3NZmz2rdvb1JSUop9PYHO0wjqgl/ydjtMmmS9Pv+93YSFvUt4+Iekpx8nJORG8vL+TmTkndxwQxhz5lh7eQpCoIFBqfJIRFYaY9r7tG9JgkZ5o0HDM28jqF9+2eqblE8ERo605l3q2ROWLl1FZuZbWENq8rj77j/y5JNPcfz4jW4DQH5w0uCgVMVQlKChI8KDRFFHULdpY1iwYD7wGpmZPxIRcQm9eg3hzTeHcMUVzc7u6y7XoHkIpQKXBo0A4ymp7Wu31ipVcmnWbBojR75GSkoK9erV4/XXX2fQoEFUrVrVPxellCo3NGgEEG9J7cK6tc6YcYYvv5xAcvIbrF27kSuuuIIPP/yQ+++/n0qVKvnvopRS5YoGjQDiba4mT91ab7nlNB988AmvvfYau3fvpm3btowZ8w133XUXNk1EKKUK0KARQDw1QaWlWfkF126tV1+dya5dH9Gy5evs3buXm266iU8++YQePXog4m48plJKadCosNzlLrw9mhSs3MYtt6SzceOHPP74Gxw4cIAuXbowYcIEunTposFCKVUoDRoVkKfcxaxZnkdWnz59mg8++IBXX32VgwcP0r17d7755hs6derk78tRSlUgGjQqIE+5i7lzLxxZ3a1bFh9//BmjRo1iz549dO/enZdffpmOHTv68xKUUhVUiL8LoIrOW+4if4zEsGG5HD36Ja1aXcWjjz5KkyZN+Pnnn5k3b54GDKVUsWnQKMdycyEhAUaNsn7m5lrb83MXrvJzF8YYpk6dSuvWrXnggQeoXr06s2bNYvHixXTp0qWsL0EpFWC0eaqc8jbmwlP32UsuWcxNNz3HsmXLuPrqq5kyZQp33XWXJriVUqVGg0Y55W3MRcHuszVr/sbs2cPo0mU69evX56OPPuLBBx8kNFR/vUqp0qXNU+WUt7wFWLmLdu32sWvXIAYPbs3PP89n9OjRbNq0ib/85S8aMJRSF4V+s5RT3sZcZGZm8tZbb/Hqq6+SlZXF4MGDGTFiBLVr1/ZbeZVSwUGDRjnlLm/RoYPh1KlvuOqq59i5cyd33XUXb7zxBpdffrm/i6uUChIaNMoBTzPTuuYtIiNTmDLl7/Tvv4Q2bdrwxRdfaG8opVSZ06DhZ956Sdls0L79fr77bhiff/45tWvXZvz48Tz88MM6maBSyi80aPiZp15SCQk5bN8+jhdeeIHMzEyeffZZhg8frs+0UEr5lQYNP3PXS8rhWMygQY9z4MAabr31Vt555x1atmzpnwIqpZQL7XLrZ+eP7j4APADcTG7uMb777jtmz56tAUMpVW5oTaOMeEp2x8ZChw65LF78AVlZw4EMmjQZxpo1w7n00shCz6uUUmVJg0YZ8Jbs/vXXNTgcA8nKWs7ll3fnmWfe5a9/vRLNcyulyiMNGmXAXbI7KSmTfv1G8f33/6Z69epMmDCB/v376zxRSqlyTYNGGbgw2f0T6emD+PbbLTz44IO8+eab1KxZ01/FU0opn2kivAycS3YfBR4EuiMijB79E5999pkGDKVUhaFBowzExkLz5tMRuQaYSFjYP+nceQ1Dh3b1d9GUUqpItHmqlBXsJdWhwxGeeuoJ1qz5mmbN2hAbO5vY2KizvaeUUqoi0aBRigr2koqImEZe3qPk5R3hpZdeYtiwYYSHh/u7mEopVWwaNErRuV5Sh4EnOH16EiEhUYwdO4chQ9r4u3hKKVVimtMoRamp4HAkAq2BKcBI8vJWcOKEBgylVGDQoFFKMjIySEp6HIgDagPJwPPY7WFERfm1aEopVWo0aJSClJQU2rZty6xZ79G48dNERiYj0ga73Rr5HRvr7xIqpVTp0JxGMeT3kEpJyWHLlteZPPkl6tWrx08//UTnzl3PPjgpKgrtJaWUCigaNIoov4fUsmXbyMgYACylbt17SU19j1q1qgMQH28tSikVaLR5qogSE2HJkilkZEQD64CJpKdPIimpur+LppRSF50GjSLIzMxk5Mi/cfp0X+BKIA3oT3q61RyllFKBrkRBQ0RqiMg8Ednk/On2z20R6SUiG0Rks4gM9fV4EWksIg4ReaYk5Syu3FxISIBRo+C999bToUMHkpM/ICzsWeAXoBlgzSulPaSUUsGgpDWNocBPxpgWwE/O9fOIiA0YB8QC1wD9xJqEyZfjxwCJJSxjseTnLu691/DCC5/y+OPt2bDhADNmJNKp0xvY7eGIoD2klFJBpaSJ8D5AF+frL4AFwP8V2KcDsNkYsxVARCY7j1vv7XgRuQPYChR4gnbZSEyEpKR00tMHAROBroSFfYVIA+bMQXtIKaWCUklrGnWNMfsAnD/ruNmnIbDLZX23c5vH40UkEit4vFxYAURkoIikiEjKoUOHin0hBc2bt4H09Bjga2AkMJfMzAakpVkBIj4eRoywfmrAUEoFi0JrGiLyI1DPzVvDffwMd4+iM4Uc8zIwxhjjKOxJdsaY8cB4gPbt2xd2Xp9MmTKFjz56CKgEzAW6A5q7UEqpQoOGMaa7p/dE5ICI1DfG7BOR+sBBN7vtBi5zWW8E7HW+9nR8DHC3iLwBVAPyROS0Mebdwi+p+LKzsxk6dChvvfUWHTrEEBb2LatXX3bec701d6GUCmYlzWnMAB4AXnP+nO5mn2SghYg0A/YA9wL9vR1vjLk5/2AReQlwXMyAkZsLEyfu44UX7mHHjl94/PHBvPXWf7DZwjV3oZRSLkoaNF4DvhGRR4CdQF8AEWkAfGyMiTPG5IjIYGAOYAM+Ncb86u34spSbCzfeuIqUlDiMOUVExER+/70/Ntu53IWO7lZKKYsYUyppgHKhffv2JiUlpUjHJCTAvfceJj29P1YP31bY7TBpkgYLpVRwEJGVxpj2vuwb9CPCU1MhI6MWVsK7FYCO8FZKKQ+CPmhER1tJblfaS0oppdwL+qARG2v1irLb0RHeSilViKCfGt1mQ0d4K6WUj4I+aID2klJKKV8FffOUUkop32nQUEop5TMNGkoppXymQUMppZTPNGgopZTyWUBNIyIih4Adbt6qBRwu4+KUN8F+D4L9+kHvAeg9APf3oIkxprYvBwdU0PBERFJ8nVclUAX7PQj26we9B6D3AEp+D7R5SimllM80aCillPJZsASN8f4uQDkQ7Pcg2K8f9B6A3gMo4T0IipyGUkqp0hEsNQ2llFKlQIOGUkopnwV00BCRXiKyQUQ2i8hQf5enLIjIZSLys4j8JiK/isiTzu01RGSeiGxy/qzu77JeTCJiE5FUEUlwrgfb9VcTkSki8rvz38KNQXgPnnL+H1gnIpNEpFKg3wMR+VREDorIOpdtHq9ZRIY5vx83iEhPXz4jYIOGiNiAcUAscA3QT0Su8W+pykQO8A9jzNXADcDjzuseCvxkjGkB/ORcD2RPAr+5rAfb9b8NzDbGXAW0wboXQXMPRKQh8ATQ3hjTGrAB9xL49+BzoFeBbW6v2fm9cC/Wc657Ae85vze9CtigAXQANhtjthpjsoDJQB8/l+miM8bsM8ascr4+hfVl0RDr2r9w7vYFcIdfClgGRKQRcBvwscvmYLr+S4FOwCcAxpgsY8xxgugeOIUClUUkFKgC7CXA74ExZhFwtMBmT9fcB5hsjDljjNkGbMb63vQqkINGQ2CXy/pu57agISJNgWhgOVDXGLMPrMAC1PFj0S62scBzQJ7LtmC6/ubAIeAzZxPdxyISSRDdA2PMHuBNYCewDzhhjJlLEN0DF56uuVjfkYEcNMTNtqDpXywiduA74O/GmJP+Lk9ZEZF44KAxZqW/y+JHoUBb4H1jTDSQTuA1w3jlbLfvAzQDGgCRIjLAv6Uqd4r1HRnIQWM3cJnLeiOs6mnAE5EwrIAx0Rgz1bn5gIjUd75fHzjor/JdZDcBt4vIdqwmya4iMoHguX6w/u3vNsYsd65PwQoiwXQPugPbjDGHjDHZwFSgI8F1D/J5uuZifUcGctBIBlqISDMRCcdK+Mzwc5kuOhERrLbs34wxb7m8NQN4wPn6AWB6WZetLBhjhhljGhljmmL9zucbYwYQJNcPYIzZD+wSkSudm7oB6wmie4DVLHWDiFRx/p/ohpXfC6Z7kM/TNc8A7hWRCBFpBrQAVhR2soAeES4icVjt2zbgU2PMv/xbootPRP4A/AKs5Vyb/j+x8hrfAI2x/kP1NcYUTJgFFBHpAjxjjIkXkZoE0fWLSBRWR4BwYCvwENYficF0D14G7sHqUZgK/AWwE8D3QEQmAV2wpj8/ALwITMPDNYvIcOBhrHv0d2NMYqGfEchBQymlVOkK5OYppZRSpUyDhlJKKZ9p0FBKKeUzDRpKKaV8pkFDKaWUzzRoKKWU8pkGDaWUUj77fyFUR0T9Owe0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "covar.plot(data='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "swedish-rapid",
   "metadata": {},
   "source": [
    "Save the exponential fit to a file with the default output name `Covariance_estimator.cov` for use in estimating the data covariance of downsampled data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "minimal-withdrawal",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "writing covariance output for Covariance estimator\n"
     ]
    }
   ],
   "source": [
    "covar.write2file(savedir='./')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "reduced-price",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
